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Abstract 

The  Rotating  Modulation  Collimator  (RMC)  belongs  to  a  larger  class  of  radiation 
imaging  systems  that  rely  on  either  temporal  or  spatial  modulation  of  incident  radiation  through 
collimation  to  map  the  location  of  the  incident  radiation  source.  The  strengths  of  these  detection 
systems  include  their  low  cost  and  simplicity.  A  major  drawback  is  the  collection  time  required 
for  low  radiation  intensities  due  especially  to  the  loss  of  radiation  information  resulting  from 
collimation.  One  method  of  addressing  this  drawback  for  the  RMC  is  by  applying  an  adaptive 
imaging  approach.  As  with  most  system  design  theory,  there  are  inherent  design  tradeoffs  for 
the  RMC.  By  using  different  RMC  configurations  for  the  same  radiation  environment 
observation,  these  tradeoffs  can  be  wagered  against  one  another  to  improve  overall  perfonnance. 
This  research  explores  the  effect  of  RMC  configuration  changes,  specifically  by  changing  the 
mask  design,  sampling  method,  and  the  angle  between  the  image  plane  and  the  RMC  centerline. 
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1.  (a)  An  overhead  view  of  how  an  RMC  uses  a  detector  without  position 
sensitivity  to  image  the  location  of  a  photon  source  based  on  the  rotation  of 
the  collimator.  As  the  collimator  rotates  the  masks  allow  a  different  ratio  of 
photons  incident  the  detector  to  transmit  through  the  collimator  to  the 
detector.  Provided  is  a  plot  of  the  ratio  of  transmitted  photons  vs  the 
collimator  rotation  angle  (b).  The  plotted  result  is  referred  to  as  the 

transmission  function . 8 

2.  This  shows  a  plot  of  a  transmission  function  (top)  with  seven  points  circled 
within  the  first  180  deg  of  rotation.  The  corresponding  exposed  detector 
images  are  shown  (bottom)  with  their  calculated  exposed  areas  in  cm2.  This 
shows  the  relationship  of  the  exposed  detector  face  relative  to  the  source  to 

the  number  of  counts  y  from  (2.1) . 11 

3.  This  shows  how  the  frequency  of  the  transmission  function  increases  as  the 
distance  between  the  source  and  the  RMC  centerline  increases  (a  through  c) 
and  how  the  phase  of  the  transmission  function  shifts  based  on  the  degree  of 

rotation  from  the  x  axis  of  the  image  plane  (c  through  e) . 12 

4.  This  is  an  example  of  a  reconstructed  image  using  two  identical  masks  which 
results  in  the  source  position  to  be  reconstructed  not  only  at  its  position  but 
also  at  its  mirror  or  symmetric  location  as  well.  The  centerline  of  the  RMC 
orthogonally  intersects  the  origin  of  the  image  plane  at  the  origin.  It  also 
shows  the  relationship  between  the  phase  and  frequency  of  the  transmission 

function  and  the  source  position  [12] . 13 

5.  (a)  This  is  an  example  of  an  observed  transmission  function,  which  is  the 
summation  of  each  independent  observation  at  a  specific  angle  for  each 
collimator  rotation.  The  circled  points  between  120  and  140  degrees  provide 
an  example  of  the  summation  results  for  the  corresponding  sample  angles 
below  (b).  For  this  example  the  observation  time  was  150  secs,  corresponding 
to  10  full  rotations.  The  statistical  bootstrap  routine  uses  the  multiple 
independent  observations  at  each  sample  angle  to  more  accurately  represent 

the  available  distribution  infonnation  in  the  averaged  transmission  function . 16 

6.  This  is  an  example  of  a  uniformly  sampled  signal  (1)  and  an  irregularly 

sampled  signal  (r).  The  irregularly  sampled  signal  shows  the  potential  to  have 
two  sample  locations  much  closer  together  than  the  uniformly  sampled  signal 
with  the  same  number  of  observations  which  can  help  in  identifying  higher 
frequency  components  of  the  signal  but  is  also  susceptible  to  under  sampling 
certain  regions  of  the  signal . 
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7.  These  plots  show  three  different  sampling  methods  used  with  the  same 
collection  time.  The  top  plot  shows  the  ideal  transmission  function  for  each 
sampling  method  on  a  semi  log  scale.  The  bottom  plot  shows  the  effect  of 
reducing  the  number  of  samples  for  a  fixed  collection  time  provides  greater 
distinction  between  the  peaks  and  valleys  of  the  transmission  function.  The 
method  using  360  samples  have  error  bars  that  essentially  overlap  throughout 

the  function,  resulting  in  ambiguous  transmission  functions . 22 
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the  centerline  of  the  detector.  Adjusting  the  pivot  angle  left  or  right  will  also 

slightly  change  the  distance  from  the  detector  to  the  image  plane . 23 

9.  Two  sample  maps  of  the  Fisher  Information  for  the  prototype  RMC  showing 
concentric  circles  of  intensity.  This  graphically  shows  the  motivation  behind 
developing  an  adaptive  pivoting  method  that  will  be  able  to  increase  the 
information  used  in  the  estimate  by  pivoting  the  detector  to  place  the  radiation 


source  in  a  ring  of  higher  intensity  [3] . 25 

10.  (L)  This  is  a  plot  of  the  ratio  of  the  mass  attenuation  coefficient  over  the 
density  for  lead  with  the  122  keV  photon  energy  highlighted.  (R)  This  is  a 
similar  plot  for  tungsten  [15] . 27 
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Adaptive  Imaging  Methods  using  a  Rotating  Modulation  Collimator 


I.  Introduction 

The  world  in  which  we  live  thrives  on  connections.  Be  it  via  internet,  transoceanic  flight, 
or  simple  public  transportation  these  avenues  connect  neighborhoods,  communities,  countries 
and  continents.  These  avenues  have  also  redefined  relationships  in  both  personal  and 
professional  circles.  People,  goods,  and  currency  travel  on  these  global  and  local  connections 
with  more  ease,  frequency,  and  intensity  than  at  any  time  in  history.  These  increases  in  dynamic 
exchanges  force  the  practice  of  radioactive  material  accountability  to  become  even  more  diligent 
and  effective.  In  April  2010  a  misplaced  Co-60  source  from  the  Delhi  University  Chemistry 
Department  in  India  was  found  at  a  salvage  yard  after  one  man  died  and  8  others  suffered  from 
the  effects  of  high  radiation  exposure.  While  the  source  originated  at  the  university,  school 
officials  were  not  aware  of  the  gamma  source  in  the  equipment  found  in  the  salvage  yard  [1].  As 
the  traffic  of  radiation  materials  increases,  both  authorized  and  otherwise,  the  need  for  detection 
systems  that  can  be  employed  in  a  variety  of  environments  and  with  materials  that  are  both  easily 
accessible  and  affordable  increases. 

The  problem  of  tracking  and  maintaining  strict  accountability  of  authorized  radioactive 
materials  and  detecting,  locating,  and  identifying  unauthorized  radioactive  materials  is  not  solved 
by  one  radiation  detection  system.  The  Department  of  Homeland  Security  uses  an  array  of 
systems  in  the  hopes  that  a  tiered  approach  will  have  the  versatility  and  flexibility  to  maintain 
accountability  and  thwart  adversary  efforts  to  harm  US  interests  [2],  The  Rotating  Modulation 
Collimator  (RMC)  has  been  explored  as  an  option  to  aid  in  stand-off  detection  of  radiation 


1 


sources  [3].  The  RMC  is  a  radiation  imaging  technique  first  developed  to  image  far  field  high 
energy  photons  in  space  without  the  use  of  focusing  lenses  or  a  position  sensitive  detector  [4]. 
The  intricate  principles  behind  the  RMC  technique  are  outlined  in  chapter  2.  The  “research 
RMC”  designed  by  Kowash  for  stand-off  detection  successfully  applied  the  RMC  technique  to 
image  radiation  sources  in  various  scenarios  [3].  The  research  completed  to  date  has  not  only 
shown  the  feasibility  of  the  system  but  uncovered  areas  that  can  be  further  exploited  to  improve 
not  only  the  RMC  as  an  option  to  the  stand-off  radiation  source  problem  but  also  other 
modulation  image  reconstruction  methods  as  well.  This  research  focuses  on  applying  adaptive 
imaging  techniques  to  the  research  RMC  to  improve  the  system  response  to  radiation. 
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I.  A.  Objectives 

The  research  for  this  thesis  focuses  on  applying  adaptive  imaging  techniques  to  a 
research  Rotating  Modulation  Collimator  (RMC).  These  techniques  can  be  broken  into  three 
distinct  areas  but  can  be  combined  to  achieve  maximum  perfonnance: 

1 .  Detennine  the  effect  of  alternative  RMC  transmission  function  sampling 
techniques.  The  traditional  sampling  technique  used  with  the  research  RMC  involves  a  near 
continuous  rotation  in  which  the  position  and  count  number  from  the  detector  are  recorded  at 
nominal  intervals  of  0.5  degrees  of  collimator  rotation.  This  data  is  then  sorted  into  1  degree 
bins,  resulting  in  360  unifonn  samples.  One  alternative  sampling  technique  is  to  look  at  the 
effect  of  decreasing  the  number  of  unifonn  samples  from  360  to  allow  for  a  longer  duration  at 
each  sample  angle.  Additionally,  irregular  sampling  using  fewer  than  360  samples  is  explored. 
Rather  than  uniform  spacing  between  samples,  irregular  sampling  samples  at  random  whole 
degree  locations.  The  desirable  effect  of  increasing  the  dwell  time  at  each  sample  angle  is  the 
improved  counting  statistics. 

2.  Detennine  the  effect  on  position  resolution  and  RMC  absolute  detection 
efficiency  of  changing  the  degree  of  modulation  through  the  masks  by  altering  the  mask  design 
parameters.  Traditional  masks  used  with  the  research  RMC  have  had  rectangular  slits  with  a  slit 
width  equal  to  half  the  pitch.  The  slit  shape  is  changed  to  a  trapezoid  and  a  variety  of  slit  widths 
are  explored  while  keeping  the  slat  width  constant  at  0.1  in. 

3.  Detennine  the  effect  of  pivoting  the  RMC  on  the  position  variance  through  the 
application  of  a  priori  Cramer-Rao  Lower  Bound  (CRLB)  of  the  position  variance.  The  CRLB 
provides  a  maximum  perfonnance  indicator  for  an  unbiased  estimate,  and  is  constant  for  a  fixed 
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RMC  system  configuration.  By  pivoting  the  centerline  of  the  RMC  with  respect  to  the  origin  of 
the  image  plane,  the  position  estimate  of  the  variance  should  follow  the  CRLB  of  the  position 
variance. 

I.  B.  Thesis  Overview 

This  thesis  is  divided  into  five  chapters.  The  first  chapter  has  provided  a  brief 
introduction  into  some  applications  of  the  RMC  along  with  the  research  objectives. 

The  second  chapter  discusses  the  theory  associated  with  general  RMC  operations  and  the 
theory  behind  the  image  reconstruction  used  with  the  research  RMC.  Also  included  is  the 
pertinent  theory  used  to  investigate  the  research  objectives,  mainly  sampling  theory,  mask 
design,  and  the  development  of  the  CRLB  on  the  position  estimate  variance.  Included  in  the 
second  chapter  is  a  section  on  the  Structural  Similarity  (SSIM)  index.  The  SSIM  is  a  full- 
reference  image  comparison  technique  developed  by  Zhou  Wang  used  to  measure  the  relative 
perfonnance  of  reconstructed  images  [5]. 

The  third  chapter  includes  a  description  of  the  RMC  and  processing  equipment  used  in 
the  research  along  with  discussions  relating  to  the  experimental  processes  used  for  each  of  the 
research  objectives.  Included  are  explanations  of  the  sources  used  and  calibration  processes 
required  prior  to  executing  the  experimental  processes. 

The  fourth  chapter  includes  the  calibration  results,  the  alternative  sampling  results,  the 
various  mask  design  results,  and  the  RMC  pivot  results.  Each  of  these  areas  presents  the  results 
along  with  an  analysis  of  the  pertinent  details. 

The  fifth  chapter  summarizes  the  major  results  of  the  research  and  presents  areas  in 
which  future  work  with  the  RMC  is  recommended. 
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II.  Theory  and  Background 

This  chapter  discusses  the  theory  relevant  to  the  RMC  and  the  background  for 
adaptive  techniques  of  interest  for  this  thesis.  A  brief  overview  of  the  basic  operating 
principles  behind  the  RMC  and  a  few  examples  of  past  applications  are  provided.  The 
focus  is  concentrated  on  three  areas  of  adaptive  imaging  using  the  RMC.  The  adaptive 
rotation  section  explores  different  sampling  methods  used  in  the  data  collection  and 
response  of  the  RMC.  This  discussion  begins  with  the  theory  that  explains  the  detector’s 
response  and  how  the  RMC  response  is  then  used  to  compute  an  estimate  of  the  source 
location  and  activity.  All  of  the  methods  are  discrete  but  vary  in  the  number  of  sampled 
points,  the  fixed  or  random  interval  between  sampling  locations,  or  both.  This  discussion 
includes  the  strengths  and  weaknesses  of  each  of  these  methods  and  when  their 
application  might  be  more  beneficial  to  our  RMC  application.  In  the  adaptive  pivot 
section  a  degree  of  freedom  is  introduced  allowing  for  the  entire  detection  system  to  alter 
its  angle  relative  to  the  desired  image  plane.  This  will  include  a  discussion  of  how  the 
knowledge  of  the  RMC  system  response  allows  for  pivot  adjustments  which  in  turn 
change  the  reconstructed  results.  The  adaptive  mask  section  explores  the  theory  behind 
selected  mask  designs  and  how  their  implementation  aides  in  our  adaptive  approach  when 
different  mask  combinations  are  used.  Finally,  a  method  for  reconstructed  image 
comparison  is  presented  as  a  tool  to  compare  the  results  of  reconstructed  images  from 
various  RMC  system  configurations.  This  comparison  tool  is  essential  when  attempting 
to  compare  differing  RMC  system  configurations  imaging  identical  radiation  scenes. 
When  approaching  adaptive  methods,  this  comparison  tool  allows  for  benchmarks  to  be 
created  which  provide  an  indication  of  a  needed  change  of  the  RMC  system.  Through  the 
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collection  process  as  these  benchmarks  are  reached  and  RMC  system  configuration 
changes  in  response  to  reconstructed  image,  the  collection  and  reconstruction  process 
becomes  adaptive. 

II.  A.  Rotating  Modulation  Collimator 

RMC  systems  belong  to  a  larger  class  of  radiation  imaging  systems  known  as 

mechanically  collimated  imaging  systems.  Mechanically  collimated  imaging  systems  are 
able  to  detennine  the  location  of  incoming  radiation  by  systematically  collimating  the 
radiation  through  any  number  of  means.  This  class  of  radiation  imaging  systems  includes 
pin-hole,  coded  aperture,  and  the  RMC  to  name  a  few.  While  certain  mechanically 
collimated  systems  are  designed  quite  simple,  a  critical  inherent  design  tradeoff  exists 
between  detection  efficiency  and  the  information  content  received  through  temporally 
varied  collimation  [6]. 

First  developed  by  astrophysicists  in  the  1960s  as  an  alternative  to  using  focusing 
lenses  or  position  sensitive  detectors,  Mertz  is  credited  with  the  RMC's  development  to 
image  cosmic  radiation  scenes  from  space  [7].  More  recently  the  RMC  has  been  adapted 
to  help  solve  the  orphan  source  problem  associated  with  standoff  detection  posed  by  the 
Department  of  Homeland  Security  in  2008  and  in  the  medical  imaging  field  as  part  of 
Neutron  Stimulated  Emission  Computed  Tomography  (NSECT)  [3]  [8].  The  premise  of 
a  rotating  modulation  collimator  is  that  a  non  position  sensitive  detector  is  able  to 
detennine  the  position  of  a  photon  source  by  rotating  a  pair  of  aligned  masks  that 
modulate  the  intensity  of  the  radiation  source.  Mask  design  is  discussed  in  detail  in 
II. 4.  A,  but  in  short  each  mask  is  designed  to  allow  a  fraction  of  the  incident  photons  to 
transmit  through  mask  openings  unobstructed.  If  the  masks  are  perfectly  aligned  and 
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have  identical  designs,  the  reconstructed  image  will  be  duplicated  about  a  line  of 
symmetry.  As  the  collimator  including  the  mask  pair  is  rotated,  the  fraction  of  photons 
allowed  to  transmit  through  the  collimator  to  the  detector  varies  as  shown  in  fig.  1 .  Each 
position  in  the  field  of  view  of  the  RMC  varies  the  transmission  of  photons  through  the 
collimator  uniquely.  Multiple  methods  have  been  used  to  solve  the  inverse  problem  of 
image  reconstruction  based  on  the  collected  data.  A  few  examples  are  CLEAN,  Fourier, 
and  Bessel  Functions  methods  [9]  [10].  The  method  used  in  this  research  estimates  the 
position  using  the  Maximum  Likelihood-Expectation  Maximization  (ML-EM)  method  by 
matching  the  measured  transmission  function  to  an  ideal  transmission  function,  an 
example  of  which  is  shown  in  fig  l.(b).  The  transmission  function  will  be  discussed  in 
further  detail  in  section  II.B.l  and  the  ML-EM  will  be  addressed  in  section  II. 2. B.  While 
the  RMC  technique  was  initially  applied  generically  to  photons,  this  research  will  focus 
specifically  on  gamma  ray  radiation  and  will  consider  all  sources  as  gamma  sources  for 
the  remainder  of  the  thesis. 
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Fig.  1.  (a)  An  overhead  view  of  how  an  RMC  uses  a  detector  without  position  sensitivity  to 
image  the  location  of  a  photon  source  based  on  the  rotation  of  the  collimator.  As  the 
collimator  rotates  the  masks  allow  a  different  ratio  of  photons  incident  the  detector  to 
transmit  through  the  collimator  to  the  detector.  Provided  is  a  plot  of  the  ratio  of 
transmitted  photons  vs  the  collimator  rotation  angle  (b).  The  plotted  result  is  referred 
to  as  the  transmission  function. 
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II.  B.  Adaptive  Rotation 

This  section  discusses  the  background  behind  the  adaptive  rotation  techniques 
investigated  in  this  research.  Initially  the  discussion  is  focused  on  how  the  transmission 
function  is  generated  by  the  RMC  and  how  it  is  then  used  to  reconstruct  an  image.  This 
information  is  important  for  all  RMC  operations,  as  it  provides  the  necessary  foundation 
of  RMC  image  generation.  Adaptive  rotation  refers  to  changing  the  rotation  of  the 
collimator  in  response  to  collected  data.  Because  the  discrete  sampling  of  the 
transmission  function  is  entirely  dependent  on  the  physical  rotation  of  the  collimator,  this 
section  could  just  as  well  be  titled  adaptive  sampling.  For  clarity,  the  term  adaptive 
rotation  will  be  used  to  retain  a  physical  sense  of  what  is  being  sampled.  A  discussion  of 
different  discrete  sampling  methods  of  the  transmission  function  concludes  this  section. 

II.  B.  1.  Transmission  Function  Generation 

The  transmission  function  for  the  RMC  refers  to  the  number  of  gammas 

transmitted  through  the  paired  masks  of  the  collimator  as  a  function  of  collimator  rotation 
angle.  As  the  paired  masks  rotate  together,  the  number  of  gammas  transmitted 
unperturbed  through  the  masks  changes  with  the  angle  of  rotation.  As  a  gamma  source  is 
observed  for  a  full  rotation,  the  transmission  function  formed  is  then  used  to  compute  an 
estimate  of  the  gamma  source  location  and  activity  based  on  knowledge  of  the  system 
response  contained  in  a  library  of  previously  computed  noiseless  transmission  functions 
for  each  possible  position  of  the  source.  It  is  widely  known  that  the  decay  of  radioactive 
materials  can  be  modeled  accurately  as  a  Poisson  stochastic  process  in  which  each  decay 
event  occurs  continuously  and  independently  from  one  another  [11],  The  equation  for  the 
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number  of  counts  detected  at  a  particular  collimator  position  from  transmission  function 
of  the  prototype  RMC  detector  used  in  this  research  is 

T„=A  oC-£(E)-^-[PH+(l-Pn)-A.(E)']  +  b(E)  .  (2-1) 

This  represents  the  number  of  counts  y  at  the  nth  collimator  angular  orientation  for 
a  dwell  time  of  x  and  is  dependent  on  the  energy  dependent  detector  efficiency  s(E), 
source  activity  a,  the  solid  angle  from  the  source  incident  the  detector  Q,  the  probability 
of  the  a  gamma  ray  transmitting  through  the  masks  at  the  nth  angle  Pn,  the  attenuation 
factor  X,  and  the  energy  dependent  background  effect  b(E)  [3].  According  to  (2.1),  for  a 
collection  with  a  constant  background  and  source  activity  and  a  stationary  source,  the 
only  change  from  one  angle  of  rotation  to  another  is  the  Pn  tenn  which  is  dependent  on 
the  mask  design.  More  specifically,  Pn  is  dependent  on  the  area  of  the  detector  face  that 
is  exposed  to  the  source.  Figure  2  depicts  this  relationship  between  the  exposed  detector 
face  to  the  incident  radiation  source  and  the  number  of  counts  received  at  an  rih  angle  of 
rotation.  The  seven  points  selected  are  the  local  minimums  and  maximums  through  the 
first  180  degrees  of  rotation.  The  images  below  the  transmission  function  provide  a  sense 
of  how  the  masks’  geometry  determines  the  exposed  detector  area.  It  is  important  to 
note,  the  masks  are  assumed  to  block  all  radiation  incident  to  the  solid  portions  of  the 
masks.  This  is  a  valid  assumption  based  on  the  low  energy  gammas  (122  keV)  and  the 
mask  material  attenuation  properties  used  in  this  research  presented  in  section  II.4.A. 
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Transmission  Function 
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Fig.  2.  This  shows  a  plot  of  a  transmission  function  (top)  with  seven  points  circled  within  the 
first  180  deg  of  rotation.  The  corresponding  exposed  detector  images  are  shown 
(bottom)  with  their  calculated  exposed  areas  in  cm2.  This  shows  the  relationship  of  the 
exposed  detector  face  relative  to  the  source  to  the  number  of  counts  y  from  (2.1). 


Consider  a  three  dimensional  space  that  includes  the  detector  and  a  gamma 
source,  where  the  position  of  the  source  is  (x, y,z)  with  z  as  the  known  distance  between 
the  detector  and  the  source  as  depicted  in  fig.  3.  The  source  location  in  the  x-y  image 
plane  is  unknown.  A  transmission  function  is  collected  by  the  RMC  after  an  observation 
period  in  which  the  collimator  is  rotated  discretely  with  a  fixed  step  interval  and  a  fixed 
dwell  time.  Figure  3  displays  the  relationship  between  the  source  position  and  the 
frequency  and  phase  of  the  generated  transmission  function. 
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Fig.  3.  This  shows  how  the  frequency  of  the  transmission  function  increases  as  the  distance 

between  the  source  and  the  RMC  centerline  increases  (a  through  c)  and  how  the  phase 
of  the  transmission  function  shifts  based  on  the  degree  of  rotation  from  the  x  axis  of  the 
image  plane  (c  through  e). 
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Figure  3  provides  a  detailed  development  the  effect  of  the  source  location  on  the 
phase  and  frequency  of  the  transmission  function.  The  relationship  can  be  succinctly 
shown  in  fig.  4.  The  distance  away  from  the  centerline  determines  the  frequency  of  the 
transmission  function  and  the  angle  of  rotation  from  the  x  axis  determines  the  phase  of 
the  transmission  function.  It  is  recognized  that  using  a  polar  coordinate  system  may 
provide  a  less  cumbersome  discussion,  but  recognizing  that  all  computations  were 
completed  using  the  Cartesian  coordinate  system,  it  will  remain  as  the  coordinate  system 
of  choice  for  consistency  purposes. 
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Fig.  4.  This  is  an  example  of  a  reconstructed  image  using  two  identical  masks  which  results  in 
the  source  position  to  be  reconstructed  not  only  at  its  position  but  also  at  its  mirror  or 
symmetric  location  as  well.  The  centerline  of  the  RMC  orthogonally  intersects  the 
origin  of  the  image  plane  at  the  origin.  It  also  shows  the  relationship  between  the  phase 
and  frequency  of  the  transmission  function  and  the  source  position  [12]. 
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A  few  interesting  points  can  be  made  concerning  this  relationship.  In  section 
II.' 4. B,  the  field  of  view  for  the  RMC  is  provided,  which  bounds  the  limits  of  the  image 
plane.  The  maximum  frequency  expected  of  a  sampled  transmission  function  is  easily 
detennined  by  the  field  of  view  and  knowledge  of  the  system  response  for  a  given  RMC 
configuration.  This  will  be  helpful  to  remember  in  discussing  sampling  techniques  in 
section  II.2.C. 

Figure  3  and  4  provide  a  host  of  ideal  transmission  functions  with  their  respective 
source  locations  in  the  image  plane.  It  also  shows  a  few  of  the  physical  source 
parameters  manifest  in  the  associated  transmission  function.  Both  figures  show  the 
frequency  increase  as  the  source  is  moved  radially  and  the  phase  shift  as  the  source  is 
rotated  from  the  horizontal  axis.  The  effect  of  the  source  distribution  is  not  considered 
for  this  research  because  the  same  source  was  used  for  all  measurements  and  the  point 
source  approximation  is  valid  for  a  disk  with  a  1  cm  diameter  over  300  cm  away  from  a 
detector  face  with  a  7.62  cm  diameter.  Another  effect  not  shown  in  figure  3  or  4  is  the 
effect  that  background  has  on  the  transmission  function.  Provided  that  the  background  is 
constant  during  the  observation  of  the  source,  the  entire  transmission  function  will  simply 
shift  vertically  according  to  the  background  intensity.  Time  varying  background  rates  are 
beyond  the  scope  of  this  research  but  may  need  to  be  incorporated  for  applications 
outside  the  lab  environment.  The  observation  period  required  to  develop  a  transmission 
function  with  both  a  distinct  phase  and  frequency  to  compare  to  the  library  of 
transmissions  functions  is  dependent  on  the  parameters  from  (2.1). 
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II.  B.  2.  Image  Reconstruction 

Upon  collection  of  an  “observed”  transmission  function,  the  next  step  is  the 
application  of  the  statistical  bootstrap  routine.  The  bootstrapping  routine  is  an  iterative 
process  that  provides  an  “averaged”  transmission  function  from  a  number  of  created 
“pseudo”  transmission  functions.  These  “pseudo”  transmission  functions  are  generated 
based  on  the  observed  data.  It  is  important  to  note  the  “observed”  transmission  function 
is  not  the  transmission  function  that  is  used  in  the  reconstruction  process.  The 
“observed”  transmission  function  is  nothing  more  than  the  summation  of  the  independent 
transmission  functions  for  every  rotation  of  the  collimator.  This  provides  an  opportunity 
to  apply  a  statistical  bootstrapping  routine  in  an  effort  to  accurately  account  for  the 
distribution  infonnation  available  for  a  single  sample  angle.  For  example,  fig.  6  (a) 
shows  an  “observed”  transmission  function  taken  over  150  seconds  with  ten  full 
rotations.  To  simplify  the  example,  the  points  between  120  and  140  are  circled  for 
further  inspection.  Figure  6(b)  shows  these  twenty  points  with  their  associated  ten 
independent  observations,  one  for  each  full  collimator  rotation.  Due  to  background  and 
the  stochastic  properties  of  radioactive  decay,  not  every  observation  at  a  given  sample 
angle  is  equal.  To  create  a  “pseudo”  transmission  function  from  the  example  in  Fig.  6, 
ten  values,  representing  the  ten  rotations,  are  selected  with  replacement  for  each  sample 
angle.  One  "pseudo"  transmission  function  represents  a  possible  transmission  function 
based  on  the  collected  data  distribution  of  each  sample  angle.  To  accurately  account  for 
the  distribution  for  all  sample  angles,  multiple  “pseudo”  transmission  functions  are 
generated.  Once  all  the  “pseudo”  transmission  functions  are  generated,  their  average  is 
calculated,  and  the  “averaged”  transmission  function  is  used  for  further  reconstruction. 
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Observed  Transmission  Function  (150  sec) 
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Observed  Transmission  Function  Contributions  for  each  full  Rotation 
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Fig.  5.  (a)  This  is  an  example  of  an  observed  transmission  function,  which  is  the  summation  of 
each  independent  observation  at  a  specific  angle  for  each  collimator  rotation.  The 
circled  points  between  120  and  140  degrees  provide  an  example  of  the  summation 
results  for  the  corresponding  sample  angles  below  (b).  For  this  example  the 
observation  time  was  150  secs,  corresponding  to  10  full  rotations.  The  statistical 
bootstrap  routine  uses  the  multiple  independent  observations  at  each  sample  angle  to 
more  accurately  represent  the  available  distribution  information  in  the  averaged 
transmission  function. 


Having  computed  an  “average”  transmission  function,  the  next  step  is  to 
reconstruct  the  image  estimate.  Having  completed  the  bootstrap  routine,  the  data  will 
simply  be  referred  to  as  the  transmission  function.  Kowash  discusses  a  few  of  the 
benefits  of  using  different  estimation  methods  for  computing  estimates  for  source 
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position  and  activity  based  on  the  collection  of  noisy  observations  [3].  By  using  the 
Maximum  Likelihood  Estimate  (MLE)  of  the  transmission  function  given  the  observed 
data,  a  source  position  and  source  activity  can  be  estimated. 

The  MLE  is  defined  for  this  problem  as, 

&ml  =  arg  q  max  L  (y  J  6) .  (2.2) 

9 ml  represents  the  parameter  values  for  0  (activity  and  position  in  this  case)  that 
maximize  the  log  likelihood  L  of  the  transmission  function  yn  from  (2.1).  L  is  defined  as 
the  log  likelihood  function  of  the  probability  density  function  of  yn  given  0 , 


L  = 


In 


p(y„\9 


(2.3) 


In  other  words  it  answers  the  question,  “What  value  of  0  (which  contains  the 
source  position  and  activity  parameters)  results  in  the  greatest  probability  of  a 
transmission  function  being  generated  matching  the  collected  data?”  The  MLE  is  often 
used  for  obtaining  practical  results.  It  can  be  easily  applied  to  complex  estimation 
problems  and  is  asymptotically  Gaussian,  unbiased,  and  efficient  for  large  enough  data 
records  [13].  While  the  MLE  works  for  estimating  multiple  parameters,  it  becomes  either 
computationally  expensive  or  altogether  impossible  to  compute  the  likelihood  maximum 
using  a  mapping  or  first  derivative  method.  The  first  derivative  method  fails  when  a 
closed  form  solution  is  unattainable,  which  then  requires  a  mapping  of  the  log  likelihood 
function  L.  A  search  of  this  map  for  the  maximum  accounts  for  the  computationally 
expensive  drawback.  However,  an  alternative  is  available.  The  method  of  parameter 
(source  activity  and  position)  estimation  used  in  this  research  is  the  Maximum  Likelihood 
Expectation  Maximization  (ML-EM).  It  is  an  iterative  method  that  is  guaranteed  to 
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converge  albeit  to  a  local  maximum  in  some  instances  and  has  a  positivity  constraint  for 
activity  estimates.  Various  forms  of  the  ML-EM  algorithm  can  be  used  and  for  this 
research  the  following  fonn  was  used, 


— 'new  — old  '  * 


A 


■(y-i  +£)]•/«• 


(2.4) 


The  new  image  estimate  following  a  given  iteration  is  .  The  previous  image  estimate 


is  Xold  and  is  equal  to  unity  for  the  first  iteration.  The  system  matrix  is  A  ,  which  along 
with  its  transpose  A '  ,  comprise  the  known  system  response  for  a  given  RMC 
configuration.  The  transmission  function  data  and  recorded  background  activity  are  y 
and  b  respectively.  Finally  a  normalization  parameter  a  is  added  to  provide  stability  to 
the  iterative  method  [3].  The  .*  and  ./  indicate  the  multiplication  and  division  operations 
are  applied  elementally.  To  help  clarify  the  dimensions  of  the  variables  from  (2.4),  the 
equation  is  rewritten  with  the  variables  removed  and  replaced  by  their  equivalent 
dimensions  in  column  row  fonnat; 

[jV,xl]=[AT,xl].*([jV1,xJVJ].([Ar,xl]./[jV,xJV),].[Arj,xl]  +  [JV,xl]))./[jV1,xl].(2.5) 

N  is  the  total  number  of  pixels  in  the  image  plane  and  Ns  is  the  total  number  of  samples 

of  the  transmission  function.  The  .*  and  ./  indicate  the  multiplication  and  division 
operations  are  applied  elementally. 

There  are  two  different  iterations  set  by  the  user  during  image  reconstruction;  one 
for  the  number  of  bootstrap  routines  to  complete  and  one  for  the  number  of  ML-EM 
iterations  to  complete.  While  the  number  of  iterations  selected  certainly  will  change  the 
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results  of  the  reconstructed  image,  results  can  be  compared  across  the  board  by  fixing 
these  iteration  parameters. 

II.  B.  3.  Sampling  Methods 

There  are  various  ways  to  sample  the  transmission  function  by  altering  the 
rotational  profile  of  the  RMC.  Sampling  theory  is  a  complex  and  deep  field  of  study, 
however,  a  few  simple  techniques  have  been  explored  for  this  research  to  determine 
whether  more  efficient  sampling  provides  similar  quality  position  estimates  using  a 
shorter  measurement  time.  Figure  7  includes  generic  examples  of  the  sampling 
techniques  explored  on  a  continuous  signal. 


Uniform  Sampling  Irregular  Sampling 


Fig.  6.  This  is  an  example  of  a  uniformly  sampled  signal  (1)  and  an  irregularly  sampled  signal 
(r).  The  irregularly  sampled  signal  shows  the  potential  to  have  two  sample  locations 
much  closer  together  than  the  uniformly  sampled  signal  with  the  same  number  of 
observations  which  can  help  in  identifying  higher  frequency  components  of  the  signal 
but  is  also  susceptible  to  under  sampling  certain  regions  of  the  signal. 
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The  Nyquist  rate, 


/«  -2-/max,  (2.6) 
determines  the  lower  bound  on  the  sampling  rate  and  is  based  on  sampling  at  least  twice 
the  highest  frequency  component  /max  of  the  continuous  signal.  Sampling  below  the 
Nyquist  rate  can  introduce  aliasing  artifacts  that  distort  the  digitized  wavefonn  from  its 
continuous  form.  Uniform  sampling,  applies  to  the  RMC  in  two  ways.  Uniform  discrete 
sampling  involves  rotating  the  collimator  a  fixed  distance  and  dwelling  for  a  fixed 
amount  of  time.  This  is  repeated  until  an  entire  revolution  or  multiple  revolutions  are 
complete.  Unifonn  continuous  sampling  is  sampling  at  a  constant  rate,  while  the  RMC 
continues  rotation.  This  is  the  method  of  sampling  the  transmission  function  that  has 
been  applied  for  most  RMCs  presented  in  literature  [3]  [4]  [8]. 

Another  method  is  irregular  sampling  which  involves  sampling  at  random 
intervals  between  consecutive  sample  locations  as  shown  by  the  generic  example  in  fig. 

7.  Using  a  collimator  rotating  discretely  with  a  fixed  step  and  fixed  dwell  time,  radiation 
scenarios  that  are  difficult  to  detect  would  require  long  collection  times  to  statistically 
strengthen  the  sample  transmission  value  at  a  given  rotation  angle.  By  applying  the 
Whittaker-  Shannon  sampling  theorem,  which  uses  the  Nyquist  Sampling  rate  as  a 
sampling  constraint,  we  can  reduce  the  collection  time  by  sampling  at  twice  the 
maximum  frequency  of  the  signal,  which  is  equivalent  to  the  field  of  view  of  our  detector 
as  discussed  in  section  II. 4. B  [14].  The  reduction  in  the  sampling  number  has  two 
significant  benefits  that  become  apparent  with  an  example.  Figure  8  shows  the 
transmission  functions  from  a  30  second  collection  using  three  different  sampling 
methods.  The  first  two  are  uniformly  sampled  with  360  and  20  samples  respectively  with 


20 


the  same  total  collection  time.  The  third  is  irregularly  sampled  with  20  samples. 

Because  the  counts  for  each  sample  are  Poisson  distributed,  the  standard  deviation  is, 

C 7  -  yfx.  (2.7) 

Where  x  is  the  number  of  counts  recorded  at  a  sample  location.  The  lower  plot  in  fig.  8 
shows  how  the  standard  deviation  for  the  reduced  sample  sizes  greatly  increases  the 
distinction  between  the  minimums  and  maximums  of  the  transmission  function.  The 


sample  size  of  360  suffers  from  overlapping  standard  deviations  throughout  the 
transmission  function,  resulting  in  ambiguous  transmission  functions  as  the  source  is 
moved  throughout  the  image  plane. 
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Counts 


Ideal  30  sec  Transmission  Function  with  various  Sampling  Methods 


Fig.  7.  These  plots  show  three  different  sampling  methods  used  with  the  same  collection  time. 

The  top  plot  shows  the  ideal  transmission  function  for  each  sampling  method  on  a  semi 
log  scale.  The  bottom  plot  shows  the  effect  of  reducing  the  number  of  samples  for  a 
fixed  collection  time  provides  greater  distinction  between  the  peaks  and  valleys  of  the 
transmission  function.  The  method  using  360  samples  have  error  bars  that  essentially 
overlap  throughout  the  function,  resulting  in  ambiguous  transmission  functions. 
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II.  C.  Adaptive  Pivot 

This  section  discusses  an  adaptive  pivoting  method  which  takes  advantage  of  the 
RMC  system  response  information  gained  by  computing  a  mapping  of  the  Cramer  Rao 
lower  bound  on  the  variance. 

II.  C.  1.  RMC  centerline  pivot 

The  ability  to  pivot  the  detector  not  only  introduces  an  additional  degree  of 
freedom  which  allows  for  expanded  field  of  view  but  has  the  added  benefit  of  potentially 
improving  the  reconstructed  image  by  altering  the  source  location  relative  to  the 
centerline  of  the  RMC  by  altering  the  pivot  angle  of  the  RMC.  Figure  9  illustrates  the 
detector  and  source  location  relationship  as  the  RMC  pivot  angle  is  adjusted. 


Fig.  8.  This  shows  an  overhead  view  of  how  the  pivot  angle  relative  to  the  image  plane  is 

adjusted  to  place  the  radiation  source  in  a  different  location  relative  to  the  centerline  of 
the  detector.  Adjusting  the  pivot  angle  left  or  right  will  also  slightly  change  the 
distance  from  the  detector  to  the  image  plane. 
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II.  C.  2.  Cramer-Rao  lower  bound  (CRLB)  on  position  variance 

The  Cramer-Rao  Lower  Bound  (CRLB)  on  the  variance  of  the  parameter  estimate 

was  calculated  for  the  system.  The  CRLB  represents  an  ideal  potential  for  any  estimate  of 
the  parameter.  For  our  work  our  parameter  estimate  will  include  both  position  and 
activity,  which  results  in  a  covariance  matrix  that  includes  the  individual  CRLB  on  the 
variance  for  the  position  and  activity  individually.  Moreover,  it  represents  the  lowest 
variance  an  unbiased  estimator  can  achieve  and  provides  a  benchmark  for  comparison. 
While  the  total  error  is  a  combination  of  the  bias  and  variance,  the  CRLB  provides  the 
lowest  achievable  variance  without  any  bias.  This  is  not  to  say,  however,  that  a  biased 
estimate  cannot  have  a  lower  variance.  The  CRLB  for  our  problem  is  derived  as, 


CRLB  =  cov (tf)>  I  (9) 


Z  =  ln[/?(yJ0)] 


(  V#  T\ 

CRLB  =  cov(^)>  -E  —ln[p(yn\0)]  C 


(2.8) 


The  Fisher  infonnation  1(9)  is  defined  as  the  negative  expectation  of  the  second 
derivative  of  the  log  likelihood  function  L  [13].  It  is  important  to  note  that  the  CRLB  is 
specific  to  a  given  RMC  configuration. 

As  previously  discussed  the  radioactive  decay  process  can  accurately  be  modeled 
as  a  Poisson  random  process.  The  Fisher  Infonnation  is  a  quantitative  way  of  measuring 
the  amount  of  infonnation  that  an  observable  random  variable  carries  about  an  unknown 


parameter  upon  which  the  probability  of  the  random  variable  depends  [13].  For  this 
application  that  is  simply  a  measurement  of  the  infonnation  content  the  transmission 
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function  has  about  the  source  position  and  activity.  Because  the  CRLB  is  dependent  on 
the  Fisher  Information,  the  two  are  closely  related. 

The  CRLB  for  RMC  systems  has  some  unique  features  that  are  worth  discussing. 
In  Kowash’s  study  of  the  CRLB  for  various  RMC  configurations  he  found  the  Fisher 
Information  of  the  position  estimates  were  not  unifonn  within  the  field  of  view. 
Subsequently,  the  CRLB  was  also  not  unifonn.  Figure  10  shows  examples  of  the  CRLB 
of  the  position  variance  maps  for  various  RMC  configurations.  The  concentric  rings 
show  that  the  information  content  of  the  RMC  increases  when  the  source  is  located  in  a 
band  of  high  intensity  compared  to  the  regions  of  lower  intensity.  Using  this  information, 
an  initial  image  can  be  reconstructed  to  determine  if  the  RMC  is  attempting  to  locate  a 
source  in  a  relative  position  known  to  have  a  high  variance  based  on  the  CRLB  map.  By 
pivoting  the  detector  to  the  left  or  right,  to  a  relative  location  of  lower  variance  the 
position  variance  of  the  reconstructed  image  decreases. 


Fig.  9.  Two  sample  maps  of  the  Fisher  Information  for  the  prototype  RMC  showing  concentric 
circles  of  intensity.  This  graphically  shows  the  motivation  behind  developing  an 
adaptive  pivoting  method  that  will  be  able  to  increase  the  information  used  in  the 
estimate  by  pivoting  the  detector  to  place  the  radiation  source  in  a  ring  of  higher 
intensity  |3J. 
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II.  D.  Adaptive  Mask 

This  section  discusses  the  effect  of  the  mask  designs  on  the  RMC  system  response 
and  subsequent  transmission  function  generation.  The  material  and  geometric  properties 
of  these  masks  detennine  the  efficacy  of  the  RMC  at  imaging  differing  radiation 
environments.  By  understanding  the  tradeoff  between  position  resolution  and  detection 
sensitivity,  masks  designed  around  these  competing  requirements  offer  two  distinct 
capabilities. 

II.  D.  1.  Basic  Mask  Design 

From  (2.1)  p  is  the  probability  a  gamma  emitted  isotropically  from  a  radiation 
source  is  incident  on  the  top  mask  of  the  RMC  and  then  successfully  passes  through  the 
open  slits  before  arriving  at  the  front  face  of  the  detector.  The  physics  which  detennine 
whether  a  gamma  will  be  blocked  by  the  masks  or  will  interact  within  the  mask  depends 
not  only  on  the  material  and  geometric  properties  of  the  mask  but  also  the  radiation  type 
and  energy  of  the  incident  radiation.  As  a  gamma  travels  through  the  solid  portions  of 
the  mask  there  are  three  major  types  of  interaction  mechanisms  that  can  occur.  These 
interactions  are  the  photoelectric  effect,  Compton  effects,  and  pair  production.  The 
combination  of  these  interactions  leads  to  the  linear  attenuation  coefficient  //  described 
as 

JU  =  T  +  (T  +  K,  (2.9) 

which  is  a  combination  of  the  probability  per  unit  path  length  for  the  photoelectric  effect 
T ,  Compton  scatter  events  cr  ,  and  pair  production  events  K  of  gammas  being  removed 
from  the  incident  gamma  flux.  It  follows  that  the  mean  free  path  of  a  gamma  l  is  defined 
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as  the  average  distance  mono-energetic  gammas  will  travel  through  a  material  before 
encountering  an  interaction  given  by  the  equation 

X  =  ~.  (2-10) 

The  masks  used  with  the  RMC  must  be  able  to  adequately  block  the  incoming  low  energy 
gammas,  meaning  the  thickness  of  the  mask  must  be  several  times  thicker  than  the  mean 
free  path  to  account  for  the  distribution  of  actual  lengths  traveled  by  the  gammas  before 
an  interaction.  The  linear  attenuation  coefficient  for  material  is  normally  listed  in  a  more 
convenient  manner  as  a  ratio  of  the  linear  attenuation  coefficient  //  to  the  density  of  the 
material  p  to  allow  for  applications  were  non  element  materials  are  involved. 


Fig.  10.  (L)  This  is  a  plot  of  the  ratio  of  the  mass  attenuation  coefficient  over  the  density  for 

lead  with  the  122  keV  photon  energy  highlighted.  (R)  This  is  a  similar  plot  for  tungsten 

[151. 


Details  of  these  three  main  interactions  are  widely  known  and  greater  details  are 
available  through  a  number  of  sources  [11]  [16].  Figure  12  shows  that  for  the  high  z 
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material  of  the  RMC  masks  (74  or  82),  the  photoelectric  effect  will  dominate  the  gamma 
interaction  with  the  mass  for  low  energy  gammas  up  to  roughly  500  keV. 


Fig.  11.  This  figure  shows  the  dependence  of  the  cross-section  for  the  three  major  types  of 
interactions  a  gammas  ray  can  undergo  as  it  transits  mass.  This  cross  section  is 
dependent  on  the  energy  of  the  gamma  ray  and  the  material  properties  of  the  mass. 
This  research  focuses  on  low  energy  gammas  and  high  z  materials  which  is  the  region 
dominated  by  the  photoelectric  effect  [11 J. 


The  standard  mask  design  implemented  in  previous  RMC  work  is  the  Ronchi 
grating.  Wilmore  explains  that  the  Ronchi  grating  is  a  transmission  grating  with  straight, 
parallel  rulings,  the  opaque  and  transmission  portions  being  of  equal  width  [11].  A  slight 
change  in  the  Ronchi  grating  is  used  by  allowing  the  slits  (open  transmission  portion  of 
the  mask)  and  slats  (opaque  portion  of  the  mask)  to  have  different  widths.  As  the  width 
of  each  grating  decreases,  the  angular  resolution  increases  at  a  cost  of  overall  detector 
sensitivity.  This  research  explores  the  effect  of  using  masks  of  varying  grate  size  all  of 
which  are  designed  after  the  Ronchi  model.  If  the  front  and  back  masks  are  identical  in 
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material  and  geometric  design  then  the  transmission  function  will  be  symmetric  through 
360  degrees  of  rotation  and  will  provide  two  source  reconstructions.  One  at  the  source 
location  and  one  at  the  negative  complement  of  the  source  position  as  in  shown  in  any  of 
the  displayed  transmission  functions  in  fig.  2. 

Part  of  the  motivation  behind  changing  the  mask  design  was  to  increase  the 
exposed  area  of  the  detector  face  to  allow  for  improved  counting  statistics  and  absolute 
detection  efficiency  discussed  in  section  II. 4. B.  This  was  accomplished  by  changing  the 
slit  shape  from  a  rectangle  to  a  trapezoid  as  shown  in  fig.  13.  This  improvement  is 
similar  to  the  improvement  seen  from  using  the  trapezoidal  rule  in  place  of  the 
rectangular  midpoint  rule  approximation  of  a  definite  integral. 


Fig.  12.  This  is  an  example  of  two  mask  designs  with  similar  pitch.  The  pitch  is  measured  from 
the  left  edge  of  one  slit  to  the  left  edge  of  the  next  slit.  The  difference  is  in  the  top  and 
bottom  shape  of  the  slit.  The  trapezoid  shape  (r)  allows  for  more  total  slit  area 
especially  for  large  slit  widths  and  only  a  few  slits. 


II.  D.  2.  Resolution  vs  Detection  Efficiency 

The  tradeoff  between  having  very  fine  slits  for  the  gammas  to  transmit  through 
and  the  time  required  to  develop  a  transmission  function  that  enables  an  accurate  position 
estimate  force  the  constraints  to  the  novel  mask  design.  This  amplifies  an  issue  that 
RMCs  and  similar  selective  transmission  techniques  must  address.  Each  radioactive 
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decay  event  radiates  information  pertinent  to  determining  the  presence,  activity, 
geometry,  and  location  of  the  source.  By  blocking  some  of  the  incident  radiation,  the 
information  carried  by  those  radiation  events  is  not  available. 

The  resolution  of  the  RMC  is  determined  by  not  only  by  the  mask  design  but  also 
the  separation  between  the  masks.  The  equation 

©  =  ^r,  (2-11) 

2  L 

shows  the  how  the  resolution  0  is  related  to  the  pitch  width  p  and  the  mask  separation 
L  [3],  The  pitch  width  is  defined  as  the  distance  between  two  left  edges  of  consecutive 
slots  in  the  mask.  This  includes  the  slot  width  and  the  distance  between  two  slots.  This 
provides  two  degrees  of  freedom  from  an  RMC  design  perspective.  We  can  improve  the 
position  resolution  by  increasing  the  number  of  slots  on  the  masks  or  increasing  the 
separation  distance  between  the  masks. 

Resolution  improvement  comes  with  a  few  tradeoff  considerations.  As  the 
number  of  slots  increases,  we  increase  the  computational  cost  of  the  reconstruction 
process.  This  occurs  because  the  reconstruction  process  relies  on  system  responses  that 
must  be  calculated  using  numerical  geometric  modeling.  As  the  number  of  slots  on  a 
mask  increases,  the  slot  width  and  pitch  decrease  for  a  fixed  mask  radius.  With 
geometrically  smaller  and  smaller  slot  dimensions,  the  number  of  points  used  to  define  an 
entire  mask  must  also  increase  to  provide  greater  fidelity  in  our  computed  system 
response.  Computational  considerations  aside,  the  lower  bound  on  the  pitch  limit  is  the 
gamma  diffraction  limit,  which  is  well  beyond  the  scope  of  this  research  but  provides  an 
interesting  concept  for  potential  RMC  applications. 
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Another  tradeoff  that  is  necessary  to  consider  when  improving  the  position 
resolution  of  the  RMC  is  the  mask  separation.  A  competing  performance  parameter  to 
improving  the  position  resolution  by  increasing  the  mask  separation  is  its  effect  on  the 


field  of  view.  The  field  of  view  FOV  is  given  by 


(2.12) 


where  d  is  the  diameter  of  the  masks  and  L  is  the  mask  separation.  Depending  on  the 
application,  the  reduction  in  field  of  view  may  not  be  as  much  of  a  consideration 
especially  for  systems  designed  with  the  ability  to  pivot  the  RMC  to  multiple  areas  with  a 
narrow  field  of  view  compared  to  fixed  systems  that  must  rely  on  one  system  position 
with  a  wider  field  of  view.  The  panoramic  advantage  to  a  pivoting  RMC  system  is 
shown  in  fig.  14. 


Top  Down  View 


RMT 


Fig.  13.  This  depicts  the  possibility  of  taking  a  panoramic  image  by  combing  images  from 
different  pivot  angles. 
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At  the  heart  of  the  discussion  of  RMC  position  resolution  is  the  quantifiable 
difference  in  the  probability  of  gammas  leaving  the  source  and  transmitting  through 
masks  to  the  detector  (p  from  (2.1))  as  the  source  is  moved  within  the  field  of  view  of  the 

image  plane.  Small  changes  in  the  source  position  resulting  in  large  differences  in  p 
have  better  position  resolution  compared  to  the  small  changes  resulting  in  minuscule 
differences  in  p  .  These  differences  in  p  are  measured  not  only  for  a  single  RMC 
collimator  position  but  for  the  entire  transmission  function. 

While  RMC  position  resolution  is  concerned  with  differences  in  p  throughout  the 
phase  of  rotation  to  determine  one  position  from  another,  the  RMC  absolute  detection 
efficiency  is  concerned  with  the  magnitude  of  p  at  each  angular  position  or  can  be 
simplified  by  the  average  p  for  a  revolution.  The  absolute  detection  efficiency  for  any 


radiation  detector  s  ,  is  defined  as 


(2.13) 


where  n  is  the  number  of  pulses  recorded  by  the  detector  and  y  .  is  the  number 
of  gammas  emitted  by  the  source  [11].  For  RMC  applications  Nemitted  is  a  constant  value 
as  the  collimator  is  rotated  based  on  the  long  half-lives  of  the  sources  compared  to  the 
collection  time.  Knowing  that  the  source  emits  isotropically,  only  the  gammas  emitted  in 
the  solid  angle  of  the  detector  face  have  a  chance  at  being  detected.  This  is  withstanding 
effects  of  the  collimator  and  the  intrinsic  efficiency  of  the  detector.  While  the  solid  angle 
from  the  source  to  the  front  face  of  the  detector  is  constant,  the  overlap  of  the  open  area 


of  the  two  masks  which  transmit  gammas  to  the  detector  changes  with  each  new  angular 
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position  as  mentioned  in  the  previous  paragraph  in  the  Pn  discussion.  To  improve  the 

RMC  system  detection  efficiency,  the  detector  face  exposed  area  needs  to  be  maximized. 
Removing  the  collimator  achieves  this  maximum.  All  gammas  emitted  from  the  source 
in  the  direction  of  the  solid  angle  of  the  detector  travel  unimpeded  to  the  detector.  While 
this  maximizes  the  system  detection  efficiency,  it  eliminates  the  ability  to  estimate  the 
position  and  intensity  of  the  gamma  source  within  the  image  plane.  As  an  alternative  to 
removing  the  collimator,  the  slot  widths  can  be  increased  to  improve  RMC  detection 
efficiency.  Herein  lays  the  tradeoff  between  position  resolution  and  RMC  detection 
efficiency.  This  tradeoff  provides  an  opportunity  to  explore  an  adaptive  approach. 

Masks  can  be  designed  to  favor  a  high  average  Pn  at  a  cost  of  position  resolution  but  still 
capable  of  providing  some  degree  of  position  information.  Alternatively,  masks  can  be 
designed  to  favor  low  position  resolution  at  the  cost  of  system  detection  efficiency.  The 
masks  used  independently  but  with  combined  results  intend  to  benefit  from  each  of  their 
design  strengths. 

II.  E.  Image  Comparison 

This  section  discusses  the  motivation  behind  an  image  comparison  method  as  it 
applies  to  adaptive  imaging,  a  few  of  the  shortcomings  of  using  descriptive  statistics  to 
compare  images,  and  finally  the  Structural  Similarity  (SSIM)  method  that  is  used  as  a 
tool  to  compare  reconstructed  images  in  this  research. 

II.  E.  1.  Adaptive  Application 

Introducing  an  adaptive  method  into  the  data  collection  of  the  RMC  requires  a 
feedback  loop.  Essentially,  an  initial  measurement  is  collected,  the  collected  data  is 
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processed  to  reconstruct  an  image,  the  image  is  analyzed  to  detennine  if  a  system 
configuration  change  is  desired,  an  additional  measurement  is  collected,  and  the 
additional  data  is  used  to  reconstruct  an  updated  image.  The  final  image  computed  using 
the  adaptive  technique  now  needs  to  be  compared  to  images  taken  with  the  same 
collection  time  for  fixed  system  configuration  to  determine  the  adaptive  effect  on 
perfonnance.  Here  two  critical  areas  have  been  identified  as  relying  upon  an  image 
analysis  method  to  not  only  determine  if  system  changes  are  needed  during  an  adaptive 
measurement  but  also  to  assess  any  performance  change  between  a  traditional  fixed 
configuration  collection  profile  and  an  adaptive  collection  profile.  A  quantifiable  method 
for  image  quality  analysis  is  essential  to  this  research. 

II.  E.  2.  Descriptive  Statistics 

Descriptive  statistics  provide  a  quantifiable  means  of  characterizing  a  set  of 
numbers.  Common  examples  of  descriptive  statistics  are  the  mean,  median,  mode, 
standard  deviation,  variance,  range,  correlation,  and  linear  regression.  F.  J.  Anscombe 
famously  illustrates  a  significant  flaw  to  singularly  relying  upon  descriptive  statistics  to 
characterize  data  sets.  His  powerful  example  computes  a  handful  of  descriptive  statistics 
for  four  separate  data  sets.  All  four  data  sets  have  identical  means  in  the  x  and  y 
directions,  variances  in  the  x  and  y  directions,  correlations  between  x  and  y,  and  linear 
regression  lines  computed  using  double  precision.  Figure  15  is  a  plot  of  each  data  set, 
clearly  showing  the  visible  differences  in  distributions  undetectable  to  the  statistics 
mentioned  previously  [17].  The  goal  for  image  comparison  is  to  retain  the  qualitative 
value  of  descriptive  statistics  but  account  for  what  is  visibly  obvious  to  the  eye  without 
having  to  manually  analyze  every  image. 
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1  2 


Fig.  14.  These  plots  show  the  Anscombe’s  quartet  of  different  distributions  with  the  same 
values  for  a  handful  of  descriptive  statistics.  This  is  a  powerful  example  of  how 
descriptive  statistics  can  be  misleading  without  visual  inspection  of  the  data.  The  goal 
for  image  comparison  is  to  retain  the  qualitative  value  of  descriptive  statistics  but 
account  for  what  is  visibly  obvious  the  eye  [17]. 

II.  E.  3.  Structural  Similarity  Index 

In  an  attempt  to  avoid  the  pitfalls  associated  with  a  strictly  descriptive  statistical 
image  analysis  the  Structural  Similarity  (SSIM)  method  developed  by  Zhou  Wang  is 
applied  to  reconstructed  images  in  this  research  [5].  Recognizing  that  the  human  eye  is 
quite  adept  at  registering  structural  content  of  an  observed  scene,  the  SSIM  method  tried 
to  incorporate  structural  information  into  the  assessment  process.  Image  assessment 
methods  normally  fall  into  either  one  of  two  categories.  The  first  is  the  blind  or  no¬ 
reference  method  that  requires  no  information  of  the  original  distortion  free  image.  The 


35 


second  is  the  full-reference  method  which  requires  a  complete  distortion  free  reference 
image.  Applications  are  typically  concerned  with  a  high  resolution  digital  image  with  an 
associated  large  file  size  as  the  reference  image  used  to  compare  images  that  have  been 
manipulated  in  some  way  for  either  effects  or  file  size  reasons.  These  altered  images  are 
then  compared  to  the  original  reference  image  to  assess  the  retained  quality  of  the 
manipulation.  The  SSIM  method  is  a  full-reference  method  [5].  In  this  application,  in 
place  of  a  reference  image  resulting  from  data  collection  and  processing,  the  truth  data  of 
the  observation  environment  is  used  in  its  place.  All  comparisons  using  the  SSIM  will 
compare  a  processed  image  to  the  truth  image  computed  by  calculating  the  activity  and 
position  of  the  source  and  using  a  background  measurement  as  the  value  for  all  other 
pixels. 

The  SSIM  method  combines  three  image  characteristics  and  compares  these 
characteristics  in  a  combined  algorithm  to  the  reference  image.  The  carryover  from 
optical  imaging  to  radiation  imaging  require  only  a  few  definition  clarifications.  The 
luminance  of  an  optical  image  is  the  same  as  the  activity  magnitude  of  a  radiation  image 
which  is  a  combination  of  the  background  and  source  in  each  pixel.  The  contrast  of  an 
optical  image  is  defined  for  radiation  images  as  the  difference  in  activity  from  pixel  to 
pixel.  The  structural  aspects  remain  essentially  unchanged  for  optical  imaging  and 
radiation  imaging. 
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The  comparative  algorithm  for  the  SSIM  index  is, 


l{x,y) 

c(x,y)- 


2/f/f  +  C, 
Mx  +  Ml  +  C\ 
2<Tx<Ty+C2 
a;  +  ct]  +  C2 

°Xy+C 3 


(2.14) 


SSIM  =  l(x,y )  •  c(x,y)  ^  ■  s(x,y) 

This  equation  is  broken  in  the  three  components  previously  discussed  as  luminance, 
contrast  and  structure.  The  reference  data  set  is  X  which  contains  the  same  number  of 
values  as  data  set  y  the  image  being  compared  to  the  reference.  The  luminance  term 


contains  the  mean  of  X  in  jux  and  the  mean  of  y  in  //, .  The  contrast  tenn  contains  the 


standard  deviation  of  I  and  y  with  ax  and  cr,  respectively.  The  structural  tenn 
consists  of  the  standard  deviations  used  in  the  contrast  term  in  addition  to  the  covariance 
between  X  and  y  in  0"n, .  Each  of  the  three  tenns  contains  a  constant  (Cj,  C2,  and  C  0  to 

avoid  instabilities  when  the  denominator  in  each  tenn  is  close  to  zero.  The  exponential 
terms  CC  ,  f3  ,  and  y  are  parameters  used  to  adjust  the  weight  or  importance  of  each  of 
the  comparative  areas.  Their  only  constraint  is  that  they  are  positive  [5].  For 
simplification  these  terms  are  set  to  unity  for  this  research.  The  solution  to  (2.14) 
provides  an  overall  index  of  comparison  of  the  measured  or  altered  image  to  the  reference 
image  [5], 
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A  visual  example  is  helpful  to  understand  the  application.  In  fig.  16,  an  image  of 
a  ball  cap  is  distorted  in  two  different  manners.  Each  distorted  image  blacks  out  the  same 
number  of  pixels,  one  with  random  pixels  selected  and  the  other  with  random  rows 
selected.  The  change  in  luminance  is  the  same  for  both  distorted  images  while  the 
change  in  contrast  appears  to  be  greater  in  the  middle  image. 


Fig.  15.  This  is  an  example  of  how  the  Structural  Similarity  (SSIM)  index  can  be  used  to 

compare  two  distorted  images  to  an  ideal  image.  The  reference  image  (1)  is  undistorted. 
The  other  two  images  have  the  same  number  of  pixels  blacked  out  with  the  difference 
being  random  pixels  being  black  out  (m)  verse  random  rows  being  blacked  out  (r). 

From  (2.12)  with  equal  exponential  weights  and  the  constants  removed  the  computed 
indexes  are  0.9808  (m)  and  0.9904  (r). 
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III.  Experiment  Setup  and  Process 

This  chapter  explores  the  experimental  equipment  and  processes  to  meet  the  research 
objectives.  Included  are  descriptions  of  the  RMC  system,  the  source  position  mechanism,  the 
radiation  sources,  and  the  software  used  for  data  processing.  Also  included  are  the  calibration 
procedures  and  the  actual  data  collection  procedures.  Finally  the  post  processing  procedures  for 
the  three  adaptive  research  areas  are  discussed.  Two  software  programs  are  used  to  operate  the 
RMC  and  process  the  data.  National  Intruments  Lab  VIEW  8.2  is  used  to  program  all  the  motor 
control  and  collection  profiles,  while  MATLAB  version  7.10.0.449  R2010a  is  used  to  save  the 
collected  data  and  for  all  post  processing. 

III.  A.  Equipment 

The  RMC  is  at  the  center  of  this  research  has  been  used  previously  in  other  areas  of  RMC 
research.  The  main  tenants  of  the  design  remain  unchanged  from  the  work  Kowash  conducted 
when  he  designed  the  RMC  in  use.  A  detailed  explanation  of  the  design  and  design  theory  for 
the  particular  RMC  used  for  this  research  is  available  in  [3]. 

III.  A.  1.  RMC 

The  main  components  of  the  RMC  are  the  flight  tube,  the  housing,  the  position  sensing 
mechanism,  the  rotating  stepper  motor,  the  pivoting  stepper  motor,  and  the  radiation  detector. 

All  these  components  are  combined  to  produce  an  electronic  pulse  that  will  be  further  processed 
and  is  discussed  in  the  following  section.  The  flight  tube  is  made  of  aluminum  and  contains  the 
mask  pair.  The  masks  will  be  discussed  in  further  detail  at  the  end  of  this  section.  The  back 
mask  is  fixed  in  place  using  set  screws  and  rests  flush  against  the  face  of  the  detector.  The  front 
mask  is  on  a  slide  allowing  for  easy  position  adjustments  from  one  measurement  to  the  next. 


The  front  mask  is  fastened  into  place  at  the  desired  mask  separation  distance  up  to  50  cm.  One 
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end  of  the  flight  tube  is  milled  down  to  fit  snugly  into  an  aluminum  driver  tube  contained  within 
the  aluminum  housing  and  fixed  to  the  driver  tube  with  set  screws.  This  driver  tube  fits  inside  of 
the  housing  and  is  free  to  rotate  with  the  aid  of  press  fit  stainless  steel  bearing  between  the 
housing  and  the  driver  tube.  The  housing  acts  to  hold  everything  in  position  and  keeps  the  fixed 
detector  aligned  with  the  rotating  driver  tube  and  flight  tube  assembly.  The  radiation  detector 
used  is  a  3  in  by  3  in  sodium  iodide  (Nal)  scintillation  detector  using  the  more  common  English 
units  for  size  measurement. 


Fig.  16.  This  is  a  picture  of  the  RMC  used  in  the  research  with  a  few  key  components  labeled.  A  few  of 
the  changes  to  the  RMC  from  previous  research  are  included.  The  idler  pulley  shown  near  the 
rotation  stepper  motor  was  added  to  provide  improved  belt  tension  for  more  consistent 
rotation.  The  position  encoder  was  upgraded  from  16  bits  to  32  bits.  Finally  the  RMC  was 
placed  on  a  motorized  pivot  table  to  allow  for  horizontal  pivoting. 
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The  rotational  position  of  the  flight  tube  must  be  accurately  recorded  throughout 
measurements  in  order  to  form  the  transmission  function.  This  is  achieved  using  a  Renishaw 
optical  encoder  ring  mounted  on  the  driver  tube.  The  encoder  ring  has  72,000  etched  lines 
equally  spaced  around  the  ring  with  a  magnetic  reset  post  mounted  at  one  location  on  the  ring. 

An  optical  pick-up  is  mounted  in  a  fixed  position  at  the  bottom  of  the  ring  to  observe  the  etching 
as  the  flight  tube  and  driver  tube  rotate.  As  an  etching  is  observed  a  square  pulse  is  generated 
and  counted  using  a  32-bit  encoding  chip.  The  counter  resets  each  time  the  magnetic  post  passes 
the  optical  pick-up.  Appendix  B.2  includes  a  picture  of  the  32-bit  encoder  chip. 

The  tube  is  rotated  using  a  HT23-397  Applied  Motion  stepper  motor  that  is  powered  and 
controlled  using  an  Applied  Motion  Si3540  control  box.  The  motor  is  mounted  on  top  of  the 
housing  and  linked  to  the  driver  tube  using  a  half  inch  timing  belt.  The  timing  belt  is  fed  through 
an  idler  gear  to  remove  any  belt  slack  and  allow  the  belt  to  wrap  around  more  of  the  drive  tube. 
The  stepper  motor  is  capable  of  stepping  less  than  1  degree  per  step  but  can  be  programmed  to 
step  at  variable  amounts.  The  gear  ratio  between  the  stepper  motor  to  the  drive  tube  is  3 .27: 1 . 

The  RMC  pivots  using  a  Velmex  B4836TS  rotary  table  with  a  Slo-Syn  stepper  motor. 

The  rotary  table  includes  a  magnetic  reference  point  that  is  used  as  a  homing  position.  The 
reference  point  discontinues  power  to  the  motor  when  passed  and  therefore  must  be  approached 
with  the  same  momentum  to  achieve  an  accurate  starting  position.  The  rotary  table  has  an 
accuracy  of  100  arc  seconds.  The  RMC  is  fastened  to  the  rotary  table  using  a  1  in  aluminum 
mounting  bracket  that  connects  the  rotary  table  to  the  RMC  housing. 

The  masks  used  in  this  research  are  3.81  cm  radius  masks  made  from  lead  or  tungsten. 
Each  mask  pair  explored  is  of  the  same  material  and  same  geometric  design.  As  discussed  in 
chapter  2  one  of  the  objectives  of  this  research  was  to  allow  more  radiation  to  be  transmitted  to 
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the  detector  while  still  gaining  some  valuable  position  infonnation.  In  fig.  18  different  mask 
designs  provide  varying  amount  of  radiation  to  transmit  through  to  the  detector.  The  ratio  of  the 
open  surface  area  to  total  mask  surface  area  is  provided  for  each  mask  design.  Based  on  the 
geometric  constraints  associated  with  straight  slot  edges,  the  slot  width  increase  does  not 
necessarily  coincide  with  an  increase  in  the  computed  surface  area  ratio  as  seen  when  comparing 
masks  2  and  3.  Table  1  provides  additional  mask  parameters. 


Mask  1  (68%)  Mask  2  (72%)  Mask  3  (70%) 


Mask  4  (69%) 


(in) 


Mask  5  (66%)  Mask  6  (57%)  Mask  7  (44%)  Mask  8  (37%) 


(in)  (in)  (in)  (in) 


Fig.  17.  This  a  two  dimensional  plot  of  the  masks  used  in  this  research.  All  masks  are  have  a  3.81  cm 

radius  and  0.635  cm  thickness.  The  percentage  next  to  the  mask  number  represents  the  ratio  of 
open  surface  area  to  total  surface  area  for  each  mask  design. 
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Table  1. _ Mask  Design  Parameters 


Mask  Number 

Number  of  Slots 

Slot  Width 
(cm) 

Slat  Width 
(cm) 

Material 

1 

2 

3.175 

0.254 

Lead 

2 

3 

2.032 

0.254 

Lead 

3 

4 

1 .4478 

0.254 

Lead 

4 

5 

1.1176 

0.254 

Lead 

5 

6 

0.889 

0.254 

Lead 

6 

9 

0.508 

0.254 

Lead 

7 

14 

0.254 

0.254 

Lead 

8 

8 

0.4 

0.4 

Tungsten 

III.  A.  2.  Sources  and  Scintillation  Detector 

The  radiation  sources  used  in  this  research  are  all  gamma  emitters.  To  significantly 
decrease  the  effects  of  partial  attenuation  through  the  masks  low  energy  gamma  sources  are  used 
to  take  advantage  of  the  lead  and  tungsten  mask  attenuation  properties.  The  source  primarily 
used  for  evaluation  of  the  RMC  image  reconstruction  is  a  Co-57  which  emits  a  122. 1  keV  with  a 
branching  ration  of  85%  and  a  half-life  of  271.8  days  [15].  The  certified  source  activity  was 
completed  most  recently  on  20  October  2009  as  0.84  mCi.  The  activity  of  the  source  changed 
slightly  over  the  course  of  the  research  but  the  majority  of  the  measurements  were  taken  when 
the  source  activity  was  0.33  mCi.  The  source  fits  into  a  bored  hole  in  a  block  of  aluminum 
which  is  attached  to  the  source  positioning  system. 

The  source  positioning  system  is  a  Velmex  Bi-slide  which  has  two  Slo-Syn  stepper 
motors  linked  to  jack  screws  to  control  the  position  of  the  source  platform.  These  motors  are 
powed  and  controlled  through  two  Velmex  VMX  control  boxes  networked  together  to  allow  for 


control  of  three  motors.  Two  motors  are  on  the  source  positioning  system  and  one  is  attached  to 
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the  rotary  table  to  pivot  the  RMC.  The  source  position  system  is  capable  of  positioning  the 
source  with  an  accuracy  of  5  //  m.  Figure  19  shows  the  source  positioning  system.  To  attain 
maximum  accuracy  the  system  must  be  leveled  at  every  junction  to  insure  alignment.  The  range 
of  motion  for  the  positioning  system  is  50  cm  in  the  horizontal  direction  and  25  cm  in  the  vertical 
direction.  The  origin  for  the  source  platfonn  was  positioned  to  always  place  the  source  within 
the  forth  Cartesian  quadrant  of  the  image  plane  (point  (0,0)  is  in  the  upper  left  in  fig.  20). 


Fig.  18.  These  photos  show  the  relative  size  of  the  point  source  used  in  this  research  (1)  and  the  source 
platform  on  the  source  positioning  system.  While  modeled  as  a  point  source,  the  source  is 
actually  disk  shaped  with  a  front  face  of  1cm  in  diameter.  The  source  platform  is  connected 
directly  to  the  Velmex  Bi-slide  and  has  aluminum  barriers  to  prevent  the  source  from  falling  to 
the  ground. 
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Fig.  19.  These  photos  show  the  source  platform  (1)  and  the  Velmex  Bi-slide  positioning  system  (r). 


One  of  the  advantages  of  the  RMC  is  the  simplistic  design  which  requires  only  one 
position  insensitive  detector.  For  this  research  a  3x3”  Bicron  Nal  detector  is  used.  These 
detectors  have  long  been  the  cheap  work  horse  of  detection  applications.  As  one  of  the  first 
scintillation  materials  used  for  radiation  detection  the  detection  properties  are  well  characterized. 
Known  for  their  relatively  high  detection  efficiency,  reliability,  and  economy,  Nal  detectors  are 
far  inferior  to  more  modem  detection  systems  in  their  energy  resolution.  While  this  research 
does  not  include  spectroscopic  applications,  it  is  essential  to  narrow  the  output  to  those  only  at 
the  energy  level  emitted  by  the  source.  The  energy  resolution  of  the  detector  is  60%  at  122  keV, 
which  will  determine  the  window  width  used  during  pulse  processing.  The  detector  fits  into  the 
housing  flush  with  the  back  mask  of  the  flight  tube  and  held  in  place  by  a  polypropylene  sleeve 


to  prevent  marring  or  scuffing.  Coupled  to  the  detector  is  a  photomultiplier  tube  (PMT)  used  to 
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convert  the  light  signal  from  the  scintillating  Nal  into  an  electronic  pulse.  For  more  details  on 
scintillation  radiation  detectors  and  PMTs  refer  to  Knoll  [11]. 


Fig.  20.  A  photo  of  the  rear  view  of  the  RMC  showing  the  sodium  iodide  (Nal)  detector  in  the  RMC 

housing  held  in  place  by  a  polypropylene  sleeve.  The  output  from  the  photomultiplier  tube  of 
the  Nal  detector  is  fed  into  the  pre-amplifier.  The  Velmex  source  positioning  system  is  visible 
in  the  background. 

III.  A.  3.  Pulse  Processing 

The  output  from  a  Nal  detector  is  proportional  to  the  energy  deposited  in  the  detector. 

For  simplified  applications  the  RMC  is  concerned  with  the  number  of  pulses  leaving  the  detector 
equivalent  to  a  particular  energy  level.  Due  to  the  poor  resolution  of  Nal  detector  this  energy 
level  becomes  an  energy  window  around  the  energy  level  of  interest.  In  order  to  determine  how 
large  to  make  the  energy  window  an  energy  spectrum  must  be  taken.  Figure  22  shows  the  block 
diagram  of  the  data  processing.  To  acquire  an  energy  spectrum  the  block  diagram  uses  the 
multi-channel  buffer  to  send  all  pulses  received  by  the  detector  and  shaped  by  the  pre-amplifier 
and  amplifier  to  the  computer  via  the  data  acquisition  card.  The  amplifier  is  set  to  ensure  that 
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desired  energy  is  well  defined  by  the  number  of  bins  in  the  energy  spectrum.  The  pulse  height 
information  is  then  processed  using  Gamma  Vision  software  with  1024  bins  to  collect  an  energy 
spectrum.  A  collection  is  gathered  until  a  visibly  recognizable  peak  is  present.  The  Gaussian 
shaped  peak  should  be  centered  about  the  channel  number  corresponding  to  the  energy  level  of 
the  selected  source.  The  channel  numbers  corresponding  to  the  points  where  the  Gaussian  curve 
meets  the  flat  background  on  either  side  of  the  center  of  the  peak  are  used  to  detennine  the 
energy  window.  Channel  numbers  are  converted  to  pulse  heights  by  recognizing  that  Nuclear 
Instrument  Module  (NIM)  equipment  is  being  used  and  the  maximum  pulse  height  response  is  3 
volts  on  NIM  equipment.  Once  the  upper  and  lower  limits  of  the  interested  energy  window  are 
detennined,  the  single  channel  analyzer  (SCA)  is  set  accordingly. 
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Fig.  21.  This  block  diagram  traces  the  RMC  data  acquisition  process  which  includes  the  pulse 
processing  |3J. 
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ITT.  B.  Calibration 

A  thorough  calibration  of  all  the  equipment  and  experiment  area  is  necessary  to  fully 
understand  the  equipment  response  to  the  changed  parameters  of  each  RMC  experiment.  This 
includes  a  calibration  of  the  Nal  detector.  Unless  stated  otherwise,  the  source  used  for  detection 
and  imaging  is  a  0.33  mCi  Co-57  source  that  emits  a  122  keV  gamma  with  a  branching  ratio  of 
85%  and  a  136  keV  with  a  branching  ratio  of  15%.  The  source  is  place  in  the  x-y  plane  a 
distance  of  3 14  cm  away  from  the  detector. 

III.  B.  1.  Pre-Collection 

Before  calibrating  the  RMC  a  few  preliminary  calibrations  must  occur.  The  first 
calibration  required  is  the  energy  calibration  for  the  Nal  detector  response  using  the  Gamma 
Vision  software.  It  is  well  known  that  while  Nal  pulse  height  detector  response  is  proportional 
to  deposited  energy,  the  proportionality  is  not  linear  across  the  entire  energy  spectrum.  For  this 
reason  an  energy  calibration  is  essential  when  looking  at  multiple  energy  levels.  While  it  may 
not  be  absolutely  necessary  in  this  work  where  only  mono-energetic  gammas  are  being 
measured,  it  is  necessary  when  using  either  a  multi-nuclide  radiation  source  or  radiation  sources 
emitting  gammas  with  higher  energy  levels.  By  observing  a  known  source  and  recording  the 
channel  number  of  the  center  of  the  corresponding  peak,  a  calibration  file  can  be  created  in 
Gamma  Vision.  The  Co-57  source  was  used  to  identify  the  122  keV  peak. 

The  next  step  includes  the  source  positioning  equipment.  The  calibration  needed  here 
consists  mostly  of  leveling  all  of  the  travel  axes  of  the  positioning  equipment.  By  using  a  bubble 
level  at  various  points  along  the  range  of  travel,  two  axis  motion  can  be  confirmed  as  truly 
horizontal  and  vertical.  Adjustments  for  unlevel  ground  or  misaligned  cross  beams  can  be  fixed 
by  using  shims  and  readjusting  the  cross  beams  (see  fig.  20). 
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The  RMC  must  also  be  carefully  calibrated  prior  to  collection.  The  masks  are  fixed  into 
position  in  the  flight  tube  using  set  screws.  The  angular  alignment  is  modeled  as  perfect  in  phase 
with  each  other  and  this  is  achieved  by  sliding  a  rigid  metal  guide  through  corresponding  slits  of 
each  of  the  masks.  The  angle  of  the  front  mask  is  adjusted  until  the  guide  rests  evenly  on  both 
masks,  which  is  verified  using  a  level.  The  angular  uncertainty  in  the  level  is  1  degree.  Once  the 
masks  are  aligned  the  starting  position  is  recorded  by  again  placing  the  metal  guide  through  both 
masks  while  the  flight  tube  is  connected  to  the  housing  and  leveling  the  guide.  The  encoder 
position  is  read  and  the  home  position  recorded  to  ensure  all  measurements  begin  from  the  same 
position.  The  final  RMC  calibration  involves  the  rotary  table.  During  the  previous  calibration 
steps,  the  pivot  angle  of  the  RMC  can  easily  have  changed.  To  ensure  that  the  RMC  is  pointed 
directly  down  range  to  the  source  positioning  system  a  homing  routine  is  executed.  This  homing 
routine  uses  the  RMC  control  software  (Labview  v8.2)  to  pivot  the  table  to  the  reference 
position,  pivot  a  set  distance  away  from  the  reference  position,  pivot  back  to  the  reference 
position,  and  then  pivot  the  set  distance  to  point  down  range  from  the  reference  position.  The 
reference  position  is  the  magnetic  pick-up  discussed  in  the  previous  section.  The  reason  the 
RMC  pivots  to  the  reference  position  twice  is  to  account  for  different  positions  the  RMC  may  be 
pointed  when  the  home  routine  is  initiated.  To  prevent  the  changes  in  momentum  from  having 
as  large  an  effect  on  the  calibration  position,  the  reference  position  is  found  when  the  RMC  is 
pivoted  from  the  same  position  at  the  same  speed  as  the  routine  homes  on  the  reference  position 
the  second  time  [18]. 

ITT.  B.  2.  Collection 

Calibration  collection  involves  collecting  data  for  various  source  positions  to  determine 
how  well  the  reconstructed  positions  using  the  RMC  correspond  with  the  actual  positions  of  the 
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source  using  the  source  positioning  system.  The  fidelity  between  these  two  data  sets  determines 
how  well  the  RMC  is  calibrated.  If  the  RMC  centerline  extends  from  the  RMC  and  is  directly 
orthogonal  to  the  origin  of  the  image  plane  the  source  moves  in  on  the  source  positioning  system 
then  the  system  is  calibrated. 

The  calibration  is  conducted  using  the  Co-57  source  and  after  all  of  the  pre-collection 
calibrations  have  been  satisfied,  the  source  is  placed  in  18  different  locations  within  the  RMC 
field  of  view.  At  each  location  the  RMC  collects  for  10  minutes  with  mask  8  and  a  mask 
separation  of  20  cm.  The  mask  selection  was  based  on  the  familiarity  of  the  mask  8  response 
based  on  previous  work  and  the  mask  separation  ensured  high  position  resolution  while  still  able 
to  span  the  range  of  motion  of  the  source  positioning  system  [3].  After  data  collection,  each  data 
set  was  processed  through  the  reconstruction  algorithm  to  provide  an  estimate  of  the  1 8  source 
positions.  The  RMC  estimates  were  compared  against  the  actual  source  positions.  Any 
discrepancies  were  noted  and  appropriate  adjustments  were  made  to  compensate.  The  easiest 
adjustment  is  to  move  the  origin  of  the  source  positioning  system  the  appropriate  direction  if  the 
discrepancy  is  consistent  for  every  point. 

Another  method  of  calibration  involves  matching  the  simulated  transmission  function 
with  that  of  a  measured  one.  One  mask  pair  is  used  to  image  one  source  location.  A  long 
collection  is  required  to  develop  a  distinguishable  transmission  function  due  to  the  detectable 
count  rate  differences  from  one  angular  position  to  the  next.  The  longer  the  run  (or  higher  the 
activity  of  the  source)  the  more  unique  the  transmission  function  becomes  as  Poisson  decay 
process  is  dominated  by  the  mean  of  the  distribution.  When  the  transmission  functions  match, 
the  simulation  has  accounted  for  the  critical  aspects  of  the  physics,  geometry,  and  material 
properties  of  the  RMC  and  background  effects. 
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III.  B.  3.  Encoder  Position  Validation 

Another  factor  to  consider  when  assessing  the  calibration  of  the  RMC  system  concerns 
the  angular  position  information  being  processed  by  the  encoder.  The  first  step  in  ensuring  the 
encoder  is  set  up  properly  involving  slowly  rotating  the  flight  tube  while  observing  the  encoder 
status  light.  The  encoder  has  three  status  indications  lights  of  green,  yellow,  and  red  referring  to 
the  amount  of  information  the  optical  pick-up  is  receiving.  The  optical  pick-up  must  be 
delicately  positioned  in  order  to  receive  the  green  status  light  through  the  entire  revolution.  The 
exception  to  this  is  when  the  magnetic  reset  passes  the  optical  pick-up  the  light  will  trip  to  red 
indicating  the  reset  of  the  counter.  These  steps  ensure  that  the  encoding  circuitry  is  receiving 
valid  infonnation. 

Once  a  measurement  is  taken  there  is  further  position  validation  that  is  done  to  ensure  the 
consistent  position  infonnation  will  be  processed  in  the  reconstruction.  Often  times  throughout 
the  course  of  runs,  particularly  long  runs,  there  will  be  position  location  that  may  inhibit  accurate 
reconstruction.  The  causes  of  these  position  anomalies  will  be  discussed  further  in  the  next 
chapter.  To  correct  for  these  errant  position  readings  two  methods  are  applied.  Both  look  to 
indentify  outliers  from  the  angular  position  data.  The  rotation  stepper  motor  is  moving  with  a 
constant  step  rate  and  distance  according  to  a  distribution.  This  distribution  is  assumed  nonnal 
and  the  outliers  outside  of  3  standard  deviations  are  identified.  Upon  identification,  these  points 
are  either  corrected  back  to  the  best  fit  line  of  the  remaining  data  (the  mean  of  the  step  distance 
distribution)  or  removed  from  consideration  entirely  along  with  the  corresponding  counts. 

III.  C.  Experimental  Data  Collection 

Upon  calibration  completion  an  experimental  data  collection  process  was  identified  that 

provided  the  flexibility  to  use  the  same  data  for  multiple  experiments.  Recognizing  that 

transmission  function  generated  during  a  collection  is  the  result  of  the  summation  of  the 
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transmission  functions  collected  for  each  revolution  of  the  RMC,  one  long  measurement  can  be 
divided  into  many  shorter  measurements  provided  a  few  conditions  are  met.  The  long 
measurement  cannot  be  broken  into  a  shorter  time  increment  than  the  time  required  for  one 
revolution.  Additionally,  the  total  measurement  must  be  taken  for  a  discrete  number  of 
revolutions  otherwise  a  partial  revolution  will  result  in  improper  equality  for  all  increments  of 
time.  For  a  given  mask  pair  the  RMC  collected  data  with  the  source  at  a  single  point  (10,-16)  for 
10  minutes  with  a  constant  rotation  step  size  and  step  rate.  Each  measurement  using  the 
available  mask  pairs  were  taken  at  separation  distances  of  20,  28,  36,  and  44  cm.  Whenever  the 
masks  were  swapped,  the  calibration  procedures  for  aligning  the  masks  and  determining  the 
home  position  using  the  encoder  were  repeated.  The  source  location  was  chosen  at  random 
within  the  range  of  the  source  positioning  system,  understanding  the  position  sensitivity  of  the 
different  mask  designs  and  mask  separations.  These  measurements  were  used  in  all  three 
adaptive  imaging  experiments.  The  adaptive  mask  experiment  compared  results  for  the  different 
mask  designs  at  a  fixed  separation  distance  and  various  corresponding  collection  times.  The 
adaptive  sampling  experiment  compared  results  for  a  mask  with  all  position  and  counting 
information  included  to  results  including  only  the  desired  sample  locations  and  corresponding 
number  of  counts.  The  adaptive  pivot  experiment  used  the  collected  data  as  a  baseline  when  the 
RMC  centerline  was  aligned  with  the  calibrated  origin  of  the  source  positioning  system.  This 
experiment  did  require  additional  data  collection  at  various  RMC  pivot  angles. 

III.  D.  Post  Processing 

This  section  will  discuss  how  the  data  collected  from  the  RMC  is  processed  to  reconstruct 
an  estimate  of  the  source  position  and  activity.  The  general  process  is  consistently  used  but  each 
adaptive  method  requires  a  tailored  method  specific  to  the  specific  application. 
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III.  D.  1.  General  Processing  Procedures 

The  raw  data  from  every  run  contains  ascending  position  information  in  the  range  of  1 
through  72,000  for  every  revolution  of  the  RMC.  It  is  important  to  point  out  that  for  one 
revolution  there  are  typically  around  800  data  points  for  the  nominal  rotation  speed.  For  each 
position  point  there  is  a  corresponding  integer  representing  the  number  of  counts  recorded  from 
the  detector  at  the  gated  energy  level  from  the  beginning  of  the  collection.  The  number  can  only 
increase  with  time  and  is  not  a  difference  in  the  number  of  counts  from  each  triggered  position. 
Again,  this  speaks  to  the  simplicity  of  the  RMC  system.  There  are  two  data  arrays.  The  first 
being  the  position  array  which  increases  in  some  rotational  speed  dependent  increment  up  to 
72,000  over  and  over  for  each  full  rotation.  With  the  second  being  the  counts  array  which 
continuously  increases  with  each  position  event  simple  counts  the  number  of  events  falling  in  the 
selected  energy  window. 

The  raw  data  containing  the  position  and  count  arrays  are  saved  at  the  end  of  each 
collection  as  a  Matlab  “.mat”  file.  Also  at  the  end  of  a  collection,  Labview  invokes  a  Matlab 
function  to  convert  the  raw  data  into  a  manageable  format.  First  the  position  information  is 
converted  into  angular  position  with  the  home  position  as  0  degrees.  Next  the  count  array  is 
converted  into  the  number  of  counts  received  at  each  angular  position  by  taking  the  difference 
between  the  total  number  of  counts  at  each  position.  Finally,  these  arrays  are  compressed  into 
360  bins  each  representing  a  degree  of  rotation  and  saved  as  a  separate  “.mat”  file.  A  series  of 
Matlab  routines  which  have  been  created  during  previous  research  are  used  to  process  the  data 
[3].  An  explanation  is  provided  here,  while  the  code  is  included  in  Appendix  A. 

Before  processing  can  begin,  the  system  response  file  must  be  created  in  order  to  use  the 
ML-EM  algorithm  discussed  in  chapter  2.  This  system  response  model  uses  parameters  of  the 
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RMC  system,  the  source,  the  background,  and  image  plane  to  compute  a  matrix  of  the 
transmission  probability  from  (2.1)  for  every  possible  source  location  in  the  image  plane  and  for 
each  degree  of  rotation.  Once  this  system  response  matrix  is  computed  it  is  stored  in  a  library  to 
expedite  future  reconstructions  with  the  same  parameters.  Next  the  reconstruction  routine  is 
invoked,  which  reconstructs  an  image  estimate  based  on  the  collected  data  using  the  system 
response  matrix  with  (2.5)  according  to  the  desired  computational  effort.  The  number  of 
bootstrap  routines  and  iteration  tolerances  are  input  parameters  at  runtime  for  this  routine.  The 
output  is  a  map  of  the  source  position  and  activity. 

III.  D.  2.  Adaptive  Rotation  Data  Processing 

After  the  initial  data  collection  for  a  given  source  position  using  all  eight  were  acquired, 
the  adaptive  experiments  could  be  completed.  For  the  adaptive  rotation  the  objective  is  to 
understand  the  effect  of  changing  the  RMC  sampling  method.  For  each  system  configuration  the 
data  are  binned  into  the  360  angular  position  increments.  Two  sampling  methods  are  explored  as 
described  in  chapter  2.  Using  uniform  sampling,  the  number  of  sampling  locations  can  be 
reduced  to  observe  the  effect  of  decreasing  the  number  of  sample  locations.  This  reduced  data 
sample  is  achieved  by  harvesting  from  the  original  collection  of  data  only  the  data  located  at  the 
desired  sample  locations.  For  example,  the  original  data  has  360  sample  locations.  Reducing  the 
number  of  uniform  samples  to  120  results  in  harvesting  the  data  contained  in  every  third  bin  of 
the  original  data.  This  changes  the  sampling  from  1  sample  per  degree  to  1  sample  per  3 
degrees.  The  collection  time  must  be  altered  to  account  for  the  reduction  in  the  number  of 
samples.  The  collection  time  is  reduced  by  the  same  factor  that  the  number  of  sampling 
locations  were  reduced. 
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For  irregular  sampling  a  similar  process  occurs  with  the  exception  that  once  the  number 
of  sample  locations  is  identified,  a  random  set  of  sample  locations  is  selected  at  varying 
intervals.  The  collection  time  still  must  be  altered  based  on  the  ratio  of  original  sample  points  to 
the  desired  sample  points. 

ITT.  D.  3.  Adaptive  Mask  Data  Processing 

The  adaptive  mask  data  processing  is  similar  to  the  adaptive  sampling  process  in  that  the 
experimental  data  is  similar.  Ten  minute  measurements  were  collected  for  each  mask  of  the 
eight  mask  design  shown  in  fig.  18  and  with  a  mask  separation  distance  of  20  cm.  A  Matlab 
routine  was  created  to  divide  the  10  minute  measurement  into  30  second  increments  to  provide 
reconstructed  of  images  of  the  progress  for  each  mask  design.  The  incremental  data  was  then 
processed  with  100  bootstrap  and  ML-EM  iterations.  The  reconstructed  images  were  saved  and 
then  run  through  a  Matlab  routine  to  compute  the  SSIM  index  based  on  the  truth  reference 
image.  These  indices  were  then  saved. 

III.  D.  4.  Adaptive  Pivot  Data  Processing 

The  adaptive  pivot  data  processing  begins  with  the  mapping  of  the  CRLB  on  the  position 
estimate  variance.  This  map  provides  the  indication  of  whether  a  RMC  pivot  adjustment  is 
required.  An  initial  image  collection  is  acquired  for  30  seconds  with  Mask  3  with  a  separation 
distance  of  20  cm.  This  data  is  then  used  to  reconstruct  an  image  as  previously  discussed  with 
100  bootstrap  and  ML-EM  iterations.  The  initial  image  reconstruction  is  analyzed  and  compared 
to  the  CRLB  map.  If  the  initial  image  is  showing  up  in  a  region  of  high  variance,  the  RMC  is 
pivoted  to  place  the  source  relative  to  the  detector  in  a  region  of  lower  position  variance.  This  is 
done  by  observing  the  distance  and  direction  required  to  move  the  source  and  calculating  the 
pivot  angle  required  to  achieve  the  shift.  It  may  be  useful  to  refer  to  fig.  9  to  understand  the 
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relationship  between  the  stationary  source  and  the  change  in  relative  distance  between  the  RMC 
centerline  and  the  source  as  the  RMC  is  pivoted. 
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IV.  Results  and  Analysis 

This  chapter  includes  the  results  of  the  calibration  and  validation  processes  discussed  in 
section  III. 2.  Also  included  are  the  results  and  analysis  of  the  experiments  discussed  in  section 
II. 3  and  III.4. 

IV.  A.  Calibration 

The  energy  calibration  for  the  Nal  detector  not  only  provided  the  energy  window  limits 
for  the  SCA,  but  also  provided  an  opportunity  to  collect  a  background  measurement  within  the 
gated  limits.  Figure  23  shows  the  Gamma  Vision  results  for  the  energy  spectrum  of  the  Co-57 
source  with  the  122  keV  peak  clearly  present.  The  peak  is  located  between  channels  94  and  134 
with  corresponding  pulse  heights  of  0.92  V  and  1.31  V.  These  are  the  upper  and  lower  limits  of 
the  SCA  representing  the  energy  level  of  interest  when  using  the  Co-57  source.  The  limits 
provide  a  valid  and  repeatable  energy  window  for  the  Co-57  provided  the  high  voltage  power 
supply  and  the  amplifier  settings  remain  constant.  As  previously  discussed,  the  poor  resolution 
of  the  Nal  detectors  results  in  relatively  larger  energy  windows  compared  to  more  modem 
detectors  like  high  purity  germanium  detectors  (HPGe).  While  the  resolution  suffers  with  the 
Nal  detector,  provided  the  energy  corresponding  to  the  source  is  accurately  bounded,  modulation 
due  to  mask  rotation  is  sufficient  enough  to  discriminate  the  probability  of  detection  from  one 
angular  position  to  another.  For  applications  including  spectroscopic  identification  with  the 
RMC,  a  detector  with  better  energy  resolution  could  easily  be  included.  Using  the  RMC  to 
properly  identify  and  subsequently  preferentially  collect  the  identified  peak  requires  spectrum 
peak  recognition  logic  not  included  in  this  research. 
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Fig.  22.  The  energy  spectrum  collected  with  the  Nal  detector  used  with  the  RMC  with  a  Co-57  source. 
The  region  of  interest  is  the  122  keV  gamma  peak  of  Co-57  with  an  associated  resolution  of 


60%. 


The  background  count  rate  used  in  the  RMC  system  response  generation  includes  only 
the  background  counts  that  fall  within  the  gated  energy  window  of  the  spectrum  discussion  in  the 
previous  paragraph.  The  application  of  this  number  includes  a  few  typical  lab  assumptions.  By 
setting  a  fixed  count  rate,  it  is  assumed  that  the  count  rate  is  time  independent.  During  the 
collection  of  data  the  background,  while  still  being  a  Poisson  stochastic  process,  the  average 
count  rate  remains  the  same.  This  is  fairly  accurate  in  the  lab  setting  but  real  world  applications 
must  accurately  account  for  a  potentially  changing  background  for  a  given  scene  as  vehicles, 
personnel,  and  equipment  can  transit  the  fixed  image  plane.  Even  in  the  lab  the  background 
count  rate  is  subject  to  change  based  on  the  possibility  of  other  experiments  being  run  within  the 
same  lab  space  using  radioactive  sources.  Another  assumption  that  must  be  considered  is  the 
background  count  rate  is  not  only  time  independent  but  also  position  independent.  In  the  lab  we 
assume  that  wherever  the  RMC  is  pointed  the  background  will  be  the  same.  This  is  a  relatively 
safe  assumption,  particularly  because  the  image  plane  is  located  314  cm  away  from  the  detector 


and  only  the  RMC  is  only  pivoted  a  few  degrees  merely  shifting  the  40  by  40  cm  image  plane  by 
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10s  of  cm.  Applying  the  pivoting  RMC  to  homeland  security  or  other  radiation  accountability 
scenarios  the  RMC  would  be  pivoted  across  a  wide  angular  range,  moving  the  image  plane 
significantly.  Each  significant  shift  would  require  an  updated  background  estimate.  The 
background  measurement  was  obtained  by  observing  the  image  plane  without  a  source  present 
for  10  minutes.  The  total  number  of  counts  was  then  divided  by  the  observation  time  resulting  in 
a  background  count  rate  of  14.93  ±  0. 15  counts  per  second.  This  measurement  was  taken 
periodically  throughout  the  course  of  the  research  and  each  measurement  fell  within  the 
uncertainty  of  the  initial  background  count  rate  estimate.  The  background  count  rate  inputted 
into  the  code  during  data  processing  was  rounded  to  15  counts  per  second. 

Additionally  the  Nal  detector  efficiency  was  calculated  in  order  to  accurately  model  the 
response  of  the  RMC.  This  simple  detection  efficiency  is  important  because  it  is  energy 
dependent.  As  discussed  in  the  previous  chapter,  RMC  applications  attempting  to  estimate  the 
energy  level  of  the  emitted  radiation  through  peak  recognition  of  a  collected  energy  spectrum 
require  a  more  thorough  detection  efficiency  calculation  for  the  entire  energy  spectrum  of 
interest.  In  this  research,  the  detection  efficiency  for  the  Nal  detector  at  122  keV  is  all  that  is 
needed.  The  intrinsic  detection  efficiency  for  the  Nal  detector  used  was  calculated  to  be  30%. 

Upon  completion  of  the  pre-collection  calibration,  the  relative  positioning  between  the 
RMC  and  the  source  positioning  system  can  move  forward.  Accurately  knowing  the  physical 
position  of  the  source  relative  to  the  RMC  is  essential  in  assessing  the  perfonnance  of  the  RMC 
reconstructed  image,  particularly  the  position  estimate.  Figure  24  shows  the  results  of  the  final 
calibration.  Previous  measurements  were  taken  with  less  agreement  indicate  an  adjustment  to 
the  origin  of  the  source  position  system  was  required  to  properly  align  the  RMC  centerline  and 
the  source  positioning  system  origin.  Figure  24  also  shows  the  disagreement  between  the 
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reconstructed  position  estimate  and  the  actual  position  estimate  near  the  origin.  As  discussed  in 
section  II,  this  is  the  blind  spot  where  the  RMC  is  unable  to  accurately  provide  an  estimate.  This 
is  due  to  the  lack  of  modulation  when  the  radiation  is  steaming  directly  down  the  centerline  of 
the  RMC  as  the  flight  tube  rotates.  The  modulation  of  the  radiation  has  to  be  larger  than  the 
noise  of  the  collected  data  in  order  to  provide  an  accurate  estimate.  The  location  of  the  source 
selected  for  the  majority  of  the  experimental  measurements  is  at  (10,-16),  an  area  of  the  image 
plane  in  which  the  RMC  accurately  estimates  the  position. 
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Fig.  23.  This  plot  shows  the  results  of  the  position  calibration  for  the  RMC.  The  actual  source  location 
is  shown  in  red,  while  the  blue  dots  indicate  the  RMC  reconstructed  image  position. 


Source  Positioning  System  Calibration 
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IV.  B.  Fidelity 

Fidelity  refers  in  this  section  refers  to  the  degree  in  which  the  computed  RMC  system 
response  agrees  with  the  actual  RMC  system  response  from  the  collected  data.  The  higher  the 
fidelity  the  closer  the  two  agree.  Ideally  there  is  a  perfect  match  between  the  simulated  and 
experimental  response,  however  a  few  things  prevent  that  from  occurring.  The  most  obvious  is 
the  randomness  associated  with  the  radioactive  decay  process.  Two  subsequent  responses  for 
either  the  simulated  or  experimental  data  will  not  agree  entirely  due  to  the  random  nature  of  the 
process.  To  account  for  this  much  longer  observation  times  are  utilized  to  compare  the  simulated 
data  to  the  collected  to  suppress  the  effects  of  the  random  noise.  Specifically,  the  transmission 
functions  are  compared  between  an  experimental  collection  and  simulated  collection.  Provided 
all  the  material,  geometric,  and  physical  properties  are  accurately  accounted  for,  the  key 
parameters  of  the  transmission  functions  would  match.  The  key  transmission  function 
parameters  include  the  frequency,  phase,  amplitude,  and  DC  shift  from  the  x  axis. 

The  transmission  function  of  the  experimental  data  is  analyzed  first.  In  fig.  25,  two  areas 
of  the  transmission  function  exhibit  erratic  behavior.  Both  these  areas  occur  in  the  last  180 
degrees  of  rotation,  and  are  quite  unexpected.  Multiple  collections  with  multiple  masks, 
separation  parameters,  source  locations,  and  sources  retain  these  oscillations  in  the  same  position 
within  the  transmission  function.  This  eliminates  a  physics  based  explanation  for  the  anomalies. 
These  anomalies  are  therefore  an  artifact  of  the  mechanical  collection  process  either  due  to 
inconsistent  belt  tension  due  to  hysteresis  or  inaccurate  rotational  position  accountability  from 
the  decoder.  The  belt  hysteresis  is  a  plausible  concern  because  the  belt  was  taut  around  the  RMC 
gearing  in  the  same  position  for  over  a  year.  Due  to  the  fact  that  identical  mask  pairs  are  used, 
the  first  180  degrees  of  the  ideal  transmission  function  matches  the  second  180  degrees  of 
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rotation.  By  invoking  this  symmetric  property,  the  last  half  of  the  experimental  transmission 
function  is  a  repeat  of  the  first  half  because  the  anomalies  are  only  in  the  second  half,  they  are 
discarded.  Figure  26  shows  the  effect  of  replacing  to  the  second  half  of  the  transmission  function 
with  the  first  half  due  to  symmetry. 
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Fig.  24.  This  is  a  plot  of  the  collected  data  using  mask  3  with  a  20  cm  mask  separation  distance  and  the 
source  located  at  (-16,8).  Looking  at  the  raw  data  there  are  two  oscillation  regions  between  180 
and  210  degrees  and  between  260  and  290  degrees.  The  smoothed  data  is  the  result  of  applying 
a  third  degree  Savitzsky-Golay  filter  with  a  span  of  21.  The  smooth  data  enables  all  but 
eliminates  the  erratic  behavior  in  the  last  180  degrees,  providing  better  symmetry  between  the 
first  180  degrees  of  rotation  and  the  last  180  degrees  of  rotation. 
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Fig.  25.  This  shows  the  symmetry  that  is  invoked  to  handle  anomalies  found  in  the  second  half  of  all  of 
the  experimentally  collect  transmission  functions.  The  upper  plot  is  the  original  transmission 
function,  while  the  lower  plot  is  the  first  180  degrees  of  the  original  data  concatenated  twice. 
This  collection  was  for  mask  3  with  a  separation  distance  of  20  cm  for  15  hours  with  the  source 
located  at  (-10,7). 


With  the  symmetric  application  applied  to  the  experimental  transmission  function,  the 
fidelity  comparison  is  essentially  narrowed  to  comparing  the  first  180  degrees  of  each  function. 
Figure  27  displays  how  well  the  simulated  and  experimental  transmission  functions  match.  This 
good  agreement  is  misleading  however,  because  the  input  values  were  altered  to  compensate  for 
significant  differences  between  the  simulated  and  experimental  transmission  functions.  In  other 
words,  the  simulation  is  modeled  using  physical  parameters  from  the  research  RMC.  The 
discrepancy  implies  that  there  is  some  unresolved  issue  when  modeling  the  new  masks  that  will 
unfortunately  skew  the  results  when  using  the  model  as  a  basis  for  image  reconstruction.  Table  2 
includes  all  the  simulated  values  that  needed  to  change  to  account  for  differences  between  the 
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experiment  and  simulation.  It  is  important  to  note,  changing  the  input  parameters  of  the  detector 
efficiency,  source  activity,  and  the  background  led  to  identifying  the  cause  of  the  discrepancy 
between  simulation  and  experiment.  The  simulation  was  only  accounting  for  radiation  that  was 
transmitted  through  both  masks  in  the  flight  tube.  This  is  only  a  fraction  of  the  radiation  the 
detector  senses  from  the  source  when  the  source  is  located  outside  the  blind  spot  of  the  RMC 
centerline.  The  radiation  can  be  transmitted  through  the  side  of  the  flight  tube  and  through  the 
back  mask  of  the  detector  bypassing  the  front  mask  entirely.  The  further  away  from  the  origin 
the  source  is  located  the  greater  the  solid  angle  between  the  source  and  the  face  of  the  detector 
which  only  includes  the  back  mask. 


x  lo4 


Fig.  26.  This  plots  shows  the  fidelity  comparison  between  the  experimental  and  simulated  transmission 
function  results.  Table  2  includes  the  alterations  required  to  bring  the  simulated  results  to  align 
with  the  experimental.  This  collection  was  for  mask  3  with  a  separation  distance  of  20  cm  for 
15  hours  with  the  source  located  at  (-10,7). 
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Table  2.  Parameter  Changes  (corrects  simulation  to  match  experimental 


transmission  function) 


detector  efficiency 

source 

activity 

background 

% 

mCi 

cnt/sec 

actual 

0.3 

0.33 

14 

correction 

0.35 

0.4 

101 

IV.  C.  Adaptive  Rotation 

The  results  for  the  adaptive  rotation  section  are  all  based  on  simulated  information. 

Based  on  the  results  of  the  calibration,  implementing  a  reliable  data  collection  process  for  the 
reduced  sampling  locations  was  not  possible.  Instead,  various  observations  scenarios  are 
simulated  in  an  attempt  to  understand  the  effect  of  reducing  the  number  of  samples.  Figure  28 
shows  an  example  of  the  simulated  measured  data  compared  to  an  ideal  transmission  function  for 
each  of  the  sampling  techniques.  As  mentioned  in  section  II.2.C,  when  the  number  of  samples 
decreases  for  a  fixed  observation  time,  the  number  of  counts  per  sample  increases  based  on  the 
increase  in  dwell  time  at  each  sampling  location.  Figure  28  again  depicts  this  with  the  increase 
in  the  average  number  of  counts  from  the  top  plot  of  roughly  1 1  for  360  sample  locations 
compared  to  the  middle  and  lower  plots  with  an  average  near  200  counts  for  20  sample  locations. 
The  middle  and  lower  plots  show  the  difference  in  shape  of  the  transmission  function  when 
sampled  uniformly  versus  irregularly  for  the  same  number  of  sample  locations.  The  effect  on  the 
uncertainty  measurement  can  be  reviewed  by  referring  to  fig  8  in  section  II. 2. C. 
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Rotation  Angle  [deg] 

Fig.  27.  These  are  plots  of  the  same  transmission  function  sampled  three  different  ways  for  a  collection 
time  of  60  seconds  and  a  source  position  of  (16,10).  The  top  plot  shows  a  uniform  sampling  of 
360  points.  The  middle  plot  is  uniformly  sampled  at  20  locations.  The  bottom  plot  is 
irregularly  sampled  at  20  locations.  The  increase  in  dwell  time  at  a  sampling  location  when  the 
number  of  sampling  locations  is  decreased  is  apparent  when  looking  at  the  same  point  in  all 
three  plots. 


In  fig.  29,  the  plots  show  various  simulated  image  reconstructions  using  uniform 
sampling.  In  this  example  the  collection  time  was  fixed  at  30  seconds  using  Mask  8  with  a 
source  position  of  (16,0)  and  the  ML-EM  and  bootstrap  iterations  fixed  at  1000  and  100  hundred 


66 


respectively.  The  results  show  the  between  images  (d)  and  (e)  false  positives  or  ghost  images 
begin  to  appear  as  the  sample  size  decreases  from  30  to  10.  Additionally  (f)  shows  more  ghost 
images  than  (e)  but  they  are  not  intense  indicating  a  degree  of  improvement  as  the  sample  size  is 
decreased  from  10  to  8.  Images  (g)  through  (i)  provide  visual  results  of  under  sampling  of  the 
transmission  function. 
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Fig.  28.  Simulated  image  reconstructions  for  various  uniform  sampling  sizes  of  the  transmission 
function  observed  for  30  seconds  using  mask  8  with  a  330  pCi  point  source  at  (16,0).  The 
number  of  ML-EM  and  bootstrap  iterations  were  fixed  at  1000  and  100  respectively. 
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By  applying  the  SSIM  to  the  reconstructed  images  in  fig.  29,  the  images  can  be  compared 
to  identify  the  effects  of  reducing  the  sample  size.  Figure  30  shows  the  SSIM  indices  for  various 
sample  sizes  for  the  point  source  located  at  three  different  locations.  For  all  three  source 
positions  the  trend  as  the  sample  size  is  reduced  is  the  same.  Because  the  radial  distance  for  all 
three  points  from  the  RMC  centerline  is  16  cm  Nyquist  sample  size  limit  for  all  three  remains  29. 
Figure  30  shows  how  the  SSIM  index  remains  relatively  steady  until  getting  below  a  sample  size 
of  30.  Below  a  sample  size  of  30,  the  SSIM  index  for  all  three  positions  oscillates  below 
plummeting  orders  of  magnitude  lower.  It  should  be  noted  that  the  SSIM  indices  for  the  position 
(16,0)  are  negative  based  on  the  negative  covariance  used  to  compute  the  SSIM  index  and 
therefore  do  not  appear  on  the  log  plot  in  fig.  30. 
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Fig.  29.  This  plot  shows  the  SSIM  index  for  various  uniform  sample  sizes  for  a  330  pCi  point  source 

using  mask  8  observed  at  three  different  locations  for  30  seconds.  The  number  of  ML-EM  and 
bootstrap  iterations  were  fixed  at  1000  and  100  respectively. 


Figure  3 1  shows  the  Nyquist  calculation  for  the  transmission  function  for  a  source 

located  at  (16,  0).  By  taking  the  fast  Fourier  transfonn  of  the  transmission,  the  highest 

meaningful  frequency  can  be  identified.  Highest  meaningful  in  this  case  refers  to  frequency 

components  that  contribute  to  the  general  overall  shape  of  the  transmission  function  not  high 

frequency  noise  components.  As  shown  in  fig.  30  the  minimum  sample  size  to  satisfy  the 

Nyquist  sampling  rate  for  the  ideal  transmission  function  at  the  point  (16,0)  is  29.  This  also 

helps  to  explain  why  ghost  images  begin  to  appear  in  fig.  29  for  sample  size  less  than  30. 
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Fig.  30.  This  shows  the  calculation  for  the  minimum  number  of  samples  required  to  satisfy  the  Nyquist 
sampling  criteria  for  the  signal  of  interest.  Plot  (a)  shows  the  ideal  transmission  function  for  a 
point  source  at  (16,0)  recorded  for  30  seconds.  Plot  (b)  is  the  fast  Fourier  transform  of  the 
transmission  function.  The  highest  meaningful  frequency  is  at  approximately  0.04  cycles/deg. 


Another  sampling  method  explored  is  irregular  sampling  as  described  in  fig.  28.  The 
problem  with  irregular  sampling  when  using  randomly  sampling  locations  is  the  image 
reconstruction  results  vary  considerably  from  subsequent  runs  due  to  the  different  sampling 
locations  change.  The  quality  of  the  image  reconstruction  is  sensitive  to  not  only  the  total 
number  of  samples  but  also  their  location.  Figure  32  illustrates  this  problem  by  comparing 
multiple  runs  for  various  sample  sizes  using  both  unifonn  and  randomly  selected  irregular 
sampling  methods.  As  the  sample  size  decreases  and  the  number  of  possible  randomly  selected 
sample  locations  increases,  the  variation  from  one  run  to  next  using  randomly  selected  irregular 


70 


sampling  increases.  Meanwhile,  because  the  sample  locations  using  uniform  sampling  are  fixed 
for  a  given  sample  size,  the  variation  between  one  run  to  the  next  is  not  an  issue. 


Fig.  31.  These  SSIM  index  comparisons  illustrate  the  difficulty  in  getting  the  same  results  for 

subsequent  runs  using  irregular  sampling  relying  on  randomly  selecting  a  fixed  number  of 
collection  locations  compared  to  the  repeatability  of  uniform  sampling.  Plot  (a)  shows  the 
SSIM  indices  for  two  subsequent  uniform  sampling  runs  where  the  points  are  overlapping  for 
the  same  sample  size.  Plot  (b)  shows  the  SSIM  indices  for  three  subsequent  random  irregular 
sampling  runs  where  the  indices  vary  for  the  same  sample  size. 


Based  on  the  infonnation  provided  in  fig.  32,  it  may  be  difficult  to  understand  the 
usefulness  of  irregular  sampling.  There  is  another  was  to  apply  irregular  sampling  that  can 
provide  more  predictable  results.  Rather  than  randomly  select  the  sampling  locations  based  on 
the  desired  sampling  size,  the  sample  locations  are  selected  based  on  knowledge  of  the 
transmission  function.  This  may  seem  a  bit  backwards  because  it  means  the  source  location  is 
already  known,  otherwise  how  is  the  transmission  function  corresponding  to  the  source  position 
realized.  This  method  could  be  applied  in  looking  in  a  specific  area  within  the  field  of  view  of 
the  RMC  to  determine  the  radiation  content.  Figure  30  shows  the  manually  selected  sample 
locations  for  the  sample  sizes  of  12  from  a  transmission  function  for  a  point  source  located  at 
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(16,0).  These  points  are  selected  because  they  represent  the  local  minima  and  maxima  within  the 
transmission  function  and  capture  the  transmission  function’s  phase  and  frequency  required  to 
identify  the  source  position.  This  will  be  referred  to  as  critical  irregular  sampling  to  differentiate 
it  from  random  irregular  sampling. 


Ideal  Transmission  Function  Mask  8  30  sec  (16,0) 


Fig.  32.  This  plot  shows  the  ideal  transmission  function  along  the  12  sample  locations  carefully  selected 
to  capture  the  frequency  and  phase  information  of  the  transmission  function. 


Figure  34  shows  the  effectiveness  of  using  critical  irregular  sampling  compared  to 
uniform  sampling  for  two  different  sample  sizes.  The  performance  improvement  for  the  critical 
irregular  sampling  is  subtle  when  the  point  source  activity  is  relatively  high  as  shown  in  images 
(a)  through  (c).  When  the  source  strength  is  reduced  from  330  pCi  to  1  lpCi,  the  critical  irregular 
sampling  method  clearly  outperforms  the  both  uniform  sampling  examples  both  in  the  computed 
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SSIM  indices  and  through  visual  inspection  of  the  reconstructed  images.  This  essentially  shows 
how  critical  irregular  sampling  can  reduce  the  minimum  detectable  activity,  or  for  this 
application  minimum  activity  required  for  reconstruction  at  a  desired  threshold.  Again,  critical 
irregular  sampling  is  only  possible  when  a  known  transmission  function  is  available  or  at  least 
identified  as  the  specific  area  of  interest  within  the  field  of  view. 
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Fig.  33.  These  simulated  reconstructed  images  show  the  impact  of  applying  critical  irregular  sampling. 
For  a  source  strength  of  330  pCi  there  is  only  minor  differences  between  the  critical  irregular 
and  uniform  sampling  reconstructions  (a  through  c).  As  the  source  strength  is  decreased,  the 
performance  of  critical  irregular  sampling  is  more  apparent  as  the  SSIM  index  becomes  nearly 
four  times  larger  than  the  uniform  sample  size  of  20  and  more  than  an  order  of  magnitude 
larger  than  the  uniform  sample  size  of  360. 
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IV.  D.  Adaptive  Mask  Design 

Each  of  the  eight  mask  designs  were  used  to  image  a  source  for  ten  minutes.  By  dividing 
up  the  collected  data  into  one  minute  increments,  a  time  lapse  profile  for  the  image  formation  for 
each  mask  design  is  made.  Fig.  35  provides  an  illustrative  example  of  how  the  images  develop 
as  collection  time  increases.  This  figure  shows  only  the  position  estimate,  but  clearly  shows  how 
the  position  resolution  deteriorates  as  the  as  the  mask  design  becomes  more  and  more 
pennissive.  Another  interesting  observation  is  the  apparent  lack  of  change  using  mask  8.  A  few 
‘ghost’  pixels  are  eliminated  from  the  60  second  reconstruction  to  the  600  second,  but  overall  the 
image  remains  unchanged.  Another  interesting  observation  is  the  result  of  mask  1.  After  ten 
minutes  of  collection  with  a  0.33  mCi  source  with  a  background  of  15  counts  per  second,  mask  1 
cannot  resolve  the  frequency  component  of  the  transmission  function  indicating  the  radial 
distance  from  the  origin  as  discussed  in  chapter  2.  While  the  phase  infonnation  of  the 
transmission  function  is  correct  at  600  seconds,  mask  1  does  not  consistently  display  the  correct 
phase  information  throughout  the  incremental  development.  Many  other  observations  can  be 
made  by  visually  analyzing  the  images  in  fig.  35,  but  most  of  these  observations  are  qualitative 
and  relative.  Where  quantitative  measures  are  needed  to  construct  fair  and  consistent 
comparisons,  the  SSIM  index  discussed  in  section  II.E.3.  is  applied. 

By  transferring  each  reconstructed  image  matrix  into  an  array,  the  time  lapse  images  for 
each  mask  design  can  be  plugged  into  (2.14).  The  reference  array  contains  the  calculated 
activity  of  the  source  at  its  known  position  with  the  value  of  the  background  in  every  other 
location.  Using  identical  masks,  even  a  perfect  reconstruction  would  not  exactly  match  the  truth 
information  in  the  reference  image  because  the  symmetric  masks  reconstruct  the  source  location 
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and  its  symmetric  complement.  Figure  36  includes  the  quantitative  comparison  of  the  time  lapse 
reconstruction  for  each  of  the  eight  masks. 
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Fig.  34.  These  plots  show  the  time  lapse  image  formation  for  four  of  the  mask  designs.  The  plots  are 

organized  in  rows  by  mask  design  and  columns  by  collection  time.  The  image  plane  is  40  by  40 
cm  with  source  located  at  (10,-16).  All  mask  separations  are  20  cm  apart  and  the  reconstruction 
used  100  bootstrap  and  MLEM  iterations.  All  images  are  each  scaled  individually. 
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Fig.  35.  This  plot  shows  the  SSIM  results  of  reconstructing  data  for  time  increments  up  to  ten  minutes 
for  each  of  the  mask  designs.  From  the  SSIM  algorithm  (2.12)  the  constants  were  set  to  zero 
and  the  exponential  parameters  were  unity.  The  computation  reconstruction  effort  was  the 
same  for  every  computation  using  100  bootstrap  and  MLEM  iterations. 


Opposed  to  the  strictly  position  estimate  comparison  in  fig.  35,  fig.  36  includes  accounts 
for  differences  in  activity  estimates  as  well.  All  masks  appear  to  exponentially  approach  some 
asymptotic  limit  of  image  quality,  shown  by  the  leveling  off  of  the  SSIM  index.  The 
implications  of  this  feature  is  that  by  after  approximately  120  seconds  each  mask  design  has 
reached  its  maximum  SSIM  index  under  these  collection  conditions.  This  may  be  useful  to 
understand  that  further  collection  time  will  not  yield  more  favorable  results  without  some  change 
to  the  RMC  system.  This  indicates  a  maximum  achievable  performance  for  a  given  mask  design 
for  a  given  imaging  scenario.  An  additional  feature  that  is  of  interest  is  the  sudden  increase  in 
the  SSIM  index  for  mask  5  between  30  and  60  seconds.  This  is  by  far  the  largest  increase  for  a 
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single  mask  over  a  30  second  increase  in  collection  time.  This  may  indicate  that  the  mask  5 
design  is  somehow  optimized  for  these  collection  conditions  between  the  30  and  60  second  time 
intervals.  Additional  experiments  would  need  to  be  run  to  verify  these  feature  and  to  attempt  to 
duplicate  is  for  the  other  mask  designs  by  altering  the  collection  conditions. 

This  plot  provides  one  part  of  the  information  necessary  to  understand  the  tradeoff 
between  position  resolution  and  detection  efficiency  for  a  mask  design.  As  would  be  expected, 
the  mask  design  that  provides  the  greatest  change  in  exposed  area  of  detector  from  one  rotation 
to  the  next  provides  the  best  image  estimate  results. 

The  other  part  of  the  mask  selection  tradeoff  is  the  detection  efficiency.  The  detection 
efficiency,  as  mentioned  in  section  II.D.2,  is  a  function  of  the  number  of  counts  recorded.  The 
efficiency  does  not  change  when  the  number  of  total  counts  is  divided  by  the  collection  time  to 
provide  the  count  rate.  The  count  rates  for  each  mask,  using  the  total  number  of  counts  detected 
in  30  second  increments  for  600  seconds,  were  recorded.  While  the  count  rate  for  a  specific 
angular  position  is  different  for  each  angle,  this  total  count  rate  includes  all  the  counts  from  the 
full  rotations. 

The  mean  count  rate  for  each  of  the  masks  is  shown  in  fig.  37.  The  trend  of  the  count 
rates  follows  the  expected  theory  with  a  few  exceptions.  Referring  back  to  the  mask  surface  area 
ratio  in  section  III. A.  1,  the  ratios  did  not  steadily  decrease  with  the  mask  design  identification 
numbers.  The  average  count  rate,  however,  falls  off  steadily  from  mask  1  to  mask  8.  None  of  the 
count  rate  uncertainties  overlap  indicating  that  the  surface  area  ratio  is  not  the  ideal  indicator  of 
mask  detection  efficiency. 
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Fig.  36.  This  is  a  plot  of  the  average  count  rate  for  each  mask  design  over  a  30  second  collection  period. 


By  combining  the  information  in  fig.  36  and  fig.  37,  the  tradeoff  between  resolution  and 
detection  efficiency  is  apparent.  An  adaptive  mask  selection  approach  can  then  be  used  by 
selecting  a  mask  with  a  high  detection  efficiency  to  observe  an  unknown  scene  for  a  particular 
energy  line.  With  prior  knowledge  of  the  background,  the  RMC  will  be  able  to  detect  the 
presence  of  a  source  earlier  with  a  lower  mask  design  number.  This  initial  collection  provides 
limited  position  information,  but  does  provide  a  crude  estimate  not  provided  by  a  bare  detector. 
The  initial  position  estimate  can  then  be  refined  by  switching  the  mask  design  to  a  higher  number 
which  has  a  better  estimate  perfonnance  based  on  fig.  36. 


78 


IV.  E.  Adaptive  Pivot 

The  CRLB  of  the  source  position  variance  provides  critical  system  performance 
information.  The  CRLB  map  in  fig.  38,  indicates  vulnerabilities  of  the  RMC  system 
perfonnance  before  any  collections  are  acquired.  Obviously  regions  of  poor  perfonnance  would 
preferably  be  avoided.  An  intennediate  image  reconstruction  indicates  if  the  source  potentially 
is  located  in  a  vulnerable  region  relative  to  the  RMC  centerline.  The  most  intense  areas  of  the 
initial  reconstruction  overlap  the  areas  of  the  CRLB  cross  section  where  the  variance  is  at  a 
maximum.  By  pivoting  the  RMC  1.7  degrees  in  the  opposite  direction  of  the  image  plane  origin 
located  314  cm  away,  the  source  is  not  located  in  CRLB  variance  trough.  The  image  of  the 
pivoted  RMC  shows  the  improvement  in  relative  variance.  The  collection  times  for  each  of  the 
reconstructed  images  are  equal.  Additionally,  the  combined  results  of  the  two  images  provide  a 
better  position  estimate  than  either  of  the  configurations  with  equivalent  collection  times. 

Another  advantage  of  using  the  rotary  table  is  the  ability  to  expand  the  field  of  view. 
While  the  FOV  indeed  follows  (2. 12),  the  idea  of  taking  a  panoramic  picture  by  splicing  multiple 
narrow  field  pictures  together  can  be  applied  as  shown  previously  in  fig.  14. 
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Fig.  37.  The  initial  reconstructed  image  (upper  r)  showing  a  high  degree  of  radial  distance  variance  for 
the  source  location.  Using  the  CRLB  position  variance  map  (upper  1)  which  is  fixed  for  given 
RMC  system  configuration,  the  pivot  angle  is  slightly  adjusted  to  avoid  areas  of  high  variance. 
The  pivoted  reconstructed  image  (lower  r)  shows  the  improvement  of  position  resolution.  The 
initial  source  location  relative  to  the  RMC  centerline  is  (16,0)  and  (22,0)  after  pivoting  the 
centerline.  The  cross-section  of  the  CRLB  map  (lower  1)  shows  a  decrease  in  the  CRLB  of  the 
position  variance  of  an  order  of  magnitude  between  16  and  22. 
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V. 


Conclusions  and  Future  Work 


V.  A.  Conclusions 

The  fidelity  results  present  a  fairly  significant  problem  in  making  any  hard  and  fast 
conclusions  concerning  the  experimental  results.  Specifically,  the  disagreement  between  the 
experimental  and  simulated  transmission  functions  and  the  erratic  behavior  of  the  experimental 
transmission  function  in  the  two  areas  of  the  last  180  degrees  of  rotation  discussed  in  section 
IV. B.  While  duplicating  the  first  180  degrees  of  rotation  was  used  as  a  work  around  for  the 
erratic  contributions,  the  experimental  and  simulated  transmission  functions  remain  unmatched. 
Fortunately,  this  mismatch  is  not  in  the  phase,  frequency,  or  general  shape  of  the  transmission 
function  and  is  simply  DC  offset  from  each  other.  Essentially  a  constant  scaling  factor  is  all  that 
separates  the  experimental  and  simulation  transmission  functions,  the  erratic  behavior 
notwithstanding.  Because  the  shape  of  the  experimental  and  simulated  transmission  functions 
are  the  same,  reliable  image  reconstructions  are  possible. 

The  effect  of  decreasing  the  number  of  uniform  samples  for  a  transmission  function  was 
identified.  Provided  that  the  Nyquist  limit  does  not  exceed  the  sampling  rate,  the  number  of 
uniform  samples  can  be  decreased  from  360.  This  can  improve  the  time  to  collect  a  well 
established  transmission  function  because  the  dwell  time  for  each  sample  can  be  increased. 
Irregular  sampling  can  provide  greatly  improved  collection  time  performance  but  is  difficult  to 
implement  due  to  the  need  to  identify  the  critical  points  of  the  transmission  function.  While 
much  fewer  points  can  be  used  to  accurately  estimate  the  frequency  and  phase  of  the 
transmission  function,  the  location  of  these  critical  points  is  only  identified  once  the  source 
position  is  known.  It  was  shown  that  by  decreasing  the  number  of  sampling  locations  by  either 
uniform  sampling  or  critical  irregular  sampling  the  minimum  activity  required  to  reconstruct  an 
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image  to  a  desire  SSIM  index  threshold  can  be  reduced.  Additionally,  the  critical  irregular 
sampling  method  can  further  reduce  the  minimum  activity  even  when  compared  to  a  reduced 
uniform  sample  size.  As  discussed,  critical  irregular  sampling  can  only  be  applied  to  a  specific 
transmission  function  prior  to  collection.  While  a  purely  adaptive  approach  was  not  identified, 
the  preliminary  understanding  of  the  effect  of  reduced  uniform  sampling  and  irregular  sampling 
provides  a  starting  point  moving  forward. 

The  mask  designs  explored  in  this  research  provide  some  enlightening  results  in  the  area 
of  mask  design.  While  the  open  trapezoid  masks  have  quite  low  position  resolution,  they 
provide  general  directional  information  of  the  source.  This  can  be  quite  useful  in  applications 
needing  detection  sensitivity  over  a  wide  field  of  view.  The  tradeoff  between  position  resolution 
and  detection  efficiency  is  readily  recognizable  when  analyzing  the  SSIM  image  quality 
assessment  and  looking  at  the  count  rate  for  all  the  mask  designs.  These  results  show  a  potential 
application  in  developing  a  collection  profile  to  improve  detection  efficiency  and  position 
resolution  rather  than  a  purely  adaptive  collection  sequence.  The  feedback  loop  required  in 
adaptive  applications  requires  SSIM  index  (or  any  other  image  quality  assessment)  benchmarks 
to  be  established  for  each  mask  configuration.  During  collection,  these  images  need  to  be 
analyzed  against  the  benchmark  to  detennine  which  RMC  configuration  change  is  necessary. 

The  information  presented  here  concerning  these  mask  designs  provides  better  characterization 
of  the  different  RMC  configurations. 

While  the  bulk  of  the  research  focused  on  sampling  techniques  and  mask  design 
evaluation,  the  RMC  pivot  experiment  provided  a  proof  of  concept  for  adaptive  pivoting.  For 
the  specific  cases  explored,  pivoting  the  mask  to  a  relative  location  with  a  lower  CRLB  on  the 
position  variance  provided  improved  perfonnance  based  on  the  SSIM  index  results. 
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V.  B.  Future  Work 

There  are  a  number  of  directions  this  research  could  follow.  One  of  the  most  promising 
directions  is  in  the  area  of  sampling  theory.  Significant  improvements  in  compressed  sensing  in 
recent  years  have  shown  considerable  reduction  in  the  lower  sampling  limit,  challenging  the 
historical  Nyquist  limit.  By  drastically  reducing  the  number  of  sampled  locations  the  RMC  may 
be  able  to  overcome  the  loss  of  infonnation  due  to  collimation.  As  one  of  the  long  standing 
drawbacks  to  RMC  application  has  been  the  long  collection  times  necessary  to  image  sources 
with  activities  near  the  background.  The  transmission  function  in  the  frequency  domain  may  be 
sparse  enough  to  take  advantage  of  the  principles  behind  compressed  sensing  introduced  in  the 
following  papers  [19]  [20].  This  would  allow  for  fewer  samples  to  be  recorded. 

Following  the  characterization  of  the  effects  of  changing  the  sampling  profile,  the  mask 
design,  and  RMC  pivot  angle,  an  ideal  collection  profile  can  be  determined.  By  adding  variable 
mask  separation  distances  to  the  flexible  RMC  system  parameters,  ideal  collection  profiles  can 
be  plotted  for  various  image  scenarios.  While  these  collection  profiles  are  not  adaptive,  they 
benefit  from  a  single  image  reconstruction  process  whereas  adaptive  methods  require  multiple 
image  reconstructions  for  the  same  radiation  scene  observation. 

An  additional  direction  concerns  mask  design.  This  research  explored  masks  having  a 
constant  slit  pitch.  Rather  than  having  to  change  a  mask  configuration  with  a  different  slit  pitch, 
the  mask  design  could  include  a  variable  slit  pitch.  A  larger  slit  pitch  in  the  middle  of  the  mask 
with  smaller  slit  pitch  on  the  edges  of  the  mask  could  potentially  provide  an  ideal  compromise 
between  the  position  resolution  and  detection  efficiency.  This  would  combine  the  position 
resolution  strength  of  the  finer  slit  pitches  with  the  detector  sensitivity  of  the  larger  slit  pitches. 
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Appendices 


Appendix  A.  Matlab  Script 


A.l  Create  RMCLib  Systems.m 

%RMC  Library  Creator 

%By:  Ben  Kowash 
%Date:  21  Jun  08 
%Update :  28  Jun  08 

o, 

o 

%Description :  This  routine  Generates  RMC  system  matrices  for  later  uses  in 
%ML-EM,  Likelihood-Function,  and  fmincon  image  reconstruction  routines. 

O, 

o 

%System  matrices  'A'  contain  the  probability  that  a  source  emitted  from  a 
%pixel  in  the  source  plane  is  transmitted  through  the  RMC  to  the  detector. 

%vl . 0 :  Assumes  black  masks 

%  Assumes  two  mask  systems 

%  Naming  convention: 

%  'Sz_L  t  x  max  dx  gridRes.mat' 

O, 

o 

%v2 . 0 :  Creates  libraries  for  multiple  RMC  measurements/simulations 

%v3 . 0 :  Utilizes  Mask  Parameter  .mat  files  to  load  specific  mask  information 

%and  build  grids  according  to  slit  shape  selection 


close  all 
clear  all 
clc 

simulation  =  'yes';  %Simulation  or  Measurment 

filename  =  'Mask8-20cm' ;  %Insert  name  of  measurement  or  simulation 

%*********************F0r  Measurements  Only************* 

Home  Address  =  'F:\Lab  RMC\ ' ; 

RMC_data_path  =  [Home  Address, ;Results\24  hr  run\ ' ] ; 

RMC_data  files  =  [RMC_data_path, 'RMC  54000  -  Mask8-20cm' ] ; 
RMC_raw_f iles  =  [RMC_data_path, ' data  54000. mat']; 

RMC_home  =  8805;  %Encoder  Position  for  Home 

S^k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k  k 


%Set  RMC  System  parameters****************************** 

num_sys  =1;  %  #  of  RMCs  or  measurements  made  at  different  positions 


RMC_CL_coord  =  [0  0];  %Cartesian  coordinates  of  RMC  centerline 

num  grids  =  [4];  %  #  of  grids  .  must  be  >=2 

grid  size  =  [2E3] ;  %  Grid  resolution 

R_det  =  [3.81];  %Radius  of  detector  [cm] 

mask_path  =  [Home  Address,  RMClib\Mask  Parameters\ ' ] ; 

mask  file  =  ['Mask8']; 

mask  coord_path  =  [Home  Address ,' MaskDesign\ ', mask  file, ' .mat '] ; 
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mask^coord  =  importdata (mask  coord_path) ; 
mask  sep  =  [21.13];  %Mask  separation  [cm] 

N  samp  =  [360];  %Number  of  samples  in  the  transmission  function 

%  samp  =  sort (randi ( 360 , 1 , N  samp));  %picks  N  samp  number  of  random  sampling 
points  between  1  and  360.  "Compressed  Sensing" 

samp  =  1:360/N  samp: 360;  %picks  N  samp  at  a  fixed  interval.  N  samp  muct  be  a 
factor  of  360. 

samp  diff  var  =  var (dif f ( samp) ) ;  %if  this  is  zero  then  it  is  a  fixed  interval 
if  ~=  0  then  random  interval 


runtime=  [30*18];  % [sec] 

tau  =  [ones(l,N  samp (1, 1) ).* (runtime (1, 1) /N_samp (1, 1) )] ;  %Dwell  time  per  bin 
[sec] 

epsRMC  =  [.6];  %Detector  Intrinsic  Efficiency 


%Set 

num 

Sxyz 

[cm] 

act 

bkg 


Source  plane  grid  parameters**************************** 
src  =1; 

=  [16  0  314];  %  location  of  source  with  detector  at  (0,0,0)  (x, 

=  [330e-6*3 . 7el0] ;  %source  activity  [Bq]  Co-57  as  of  190ct2010 
=  [15] ;  %background  [cps] 


Y/ 


z) 


%Set  RMC  System  Matrix  Parameters************************** 

x  max  =  40;  %Number  of  pixels  in  x 

y_max  =  40;  %Number  of  pixels  in  y 

dx  base  =  1;  %Pixel  resolution  in  x  [cm] 

dy_base  =  dx  base;  %Pixel  resolution  in  y  [cm] 

x  expect  =  0;  %Used  to  center  window  to  zoom  in  Measurement 

y  expect  =  0;  %Used  to  center  window  to  zoom  in  Measurement 

z^max  =  1 ; 

z  0  =  314;  %location  of  reconstructed  image  plane, 
dz  =  100;  %Resolution  in  z  [cm] 


%Corrects  the  repetitive  sys  parameters  for  the  correct  num  of  sys 
if  num  sys  >  1; 

RMC_CL_coord  =  correct_f or_num_sys (num  sys, RMC  CL_coord) ; 
num  grids  =  correct  for  num^sys (num^sys , num  grids); 
grid  size  =  correct  for  num^sys (num^sys , grid  size) ; 

R_det  =  correct_for_num_sys (num^sys, R_det) ; 

mask  file  =  correct  for  num^sys (num^sys , mask  file) ; 

mask_sep  =  correct  for  num^sys (num_sys , mask  sep); 

N  samp  =  correct_f or^num  sys (num  sys,N  samp) ; 
samp  =  correct  for  num^sys (num  sys, samp); 
runtime  =  correct__for_num  sys  (num  sys,  runtime)  ; 
tau  =  correct_f or_num_sys (num_sys,  tau)  ; 
epsRMC  =  correct_f or^num  sys (num  sys, epsRMC); 

end; 


%Reads  in  mask  parameter  information  for  each  RMC  system 
slot_len ( 1 : 30 , 1 : num_sys )  =  0; 
for  i  =  1 : num  sys; 

mask  name  =  [mask_path, mask  file  (i,  :),' .mat '] ; 
load  (mask  name) ; 
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a ( i , : ) =slit_pitch;  %slit  pitch  [cm] 
w (i,  : ) =slit_width;  %slit  width  [cm] 

c (i, : ) =lef t_edge_lst_slit;  %location  of  left  edge  of  furthest  left  slit 

[cm] 

t  mask ( i ) =mask  thickness;  %mask  thickness  [cm] 

mask  min(i)=slit  to  radius;  %minimum  distance  between  slit  and  edge  of 
mask  [cm] ,  this  determines  the  slot  length  unless  specified 

slot_type (i, : ) =slit_shape;  %describes  the  shape  of  the  slit  as  string 
variable  (trapazoid,  rounded,  rectangle,  or  specified) 

if  strcmp (slot_type (i, :),' specified' ) ;  %' specified'  allows  for  specific 

slot  length  dimensions  as  input  for  a  rectangle  slit 

slot_len  ( 1 :  length  ( slot__length)  ,  i )  =slot  length; 

end; 

mu(i)=mass  atten;  %mass  attenuation  coefficent  based  on  the  material 
type  and  energy  of  gammas 
end; 


%Define  physical  properties  of  the  masks 
atten  =  exp (-mu. *mean (t  mask, 2)); 


%***  Compute  mask  separation  parameters  'L'  for  each  grid  *** 
for  i=l:num_sys 

if  (num  grids  ==  2) 

L(i,:)  =  [0,  mask^sep] ; 

else 

mask_pos  =  0; 

for  j=l:num  grids (i) 

L(i,j)  =  mask_pos; 
if  (isEven(j)  ==  0) 

mask_pos  =  mask_pos  +  t_mask(i); 

else 

mask_pos  =  mask_pos  +  mask  sep(i); 

end 

end 

end 

end 

clear  mask  pos 


%Corrects  the  repetitive  source  parameters  for  the  correct  num  of  sys 
if  num_src  >  1; 

Sxyz  =  correct  for  num  sys (num  src, Sxyz) ; 
act  =  correct_f or_num_sys  (num_src,  act)  ; 
bkg  =  correct^for  num^sys (num  src,bkg); 

end; 
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%Saves  input  into  appropriate  library  file 
if  (strcmp (simulation, ' yes ') ==1) 

sim_system  =  [Home  Address RMCLib\Sim  System\ ' ] ; 
sim_src  =  [Home  Address, ' RMCLib\Sim_Src\ '] ; 
sim_sys_name  =  [sim  system, filename, 1  simsys .mat ' ] ; 

save (sim  sys  name,  ' num  sys '  ,  ' RMC  CL  coord ' num  grids',  'grid  size '  ,  ' R  det '  ,  ' a 
' ,  ' w ' ,  ' c ' ,  ' mask  min',  ' t  mask '  ,  ' mask  sep ' ,  '  L ' ,  'slot  len '  ,  'slot  type ' ,  ' N  samp ' , 

' samp ' ,  ' tau '  ,  ' epsRMC '  ,  'mu',  ' atten '  ,  ' mask^coord ' ) 

sim_src_name  =  [sim  src, filename, '  simsrc .mat ' ] ; 
save ( sim_src_name, ' num  src ' , ' Sxyz ' , ' act ' , ' bkg ' ) 

else  %Measurement  Data 

meas  system  =  [Home  Address ,' RMCLib\Meas_System\ '] ; 
meas_sys_name  =  [meas_system, filename, '  meassys . mat ' ] ; 

save (meas  sys  name,  ' Sxyz '  ,  ' num  sys '  ,  ' RMC  CL  coord ',  ' num  grids',  'grid  size',  ' R 
det ' , ' a ' , ' w ' , ' c ' , ' mask  min ' , ' t  mask ' , ' mask  sep', ' L ' , 'slot  len', 'slot  type ' , ' 
N_samp ' , ' tau ' , ' epsRMC ' , 'mu', ' atten ' ) 

meas  data  =  [Home  Address, 'RMCLib\Meas  data\ ' ] ; 

meas_data  name  =  [meas_data, filename, '  measdata . mat ' ] ; 

save (meas  data  name, 'RMC  data  files', 'RMC  raw  files', 'RMC  home') 

end 


~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  ~k 

%*****************general-e  System  Matrix******************************** 

S-'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


%Sets  up  the  number  of  image  planes  to  look  at  seperated  by  dz 
z  =  zeros (z  max, 1) ; 
for  i=l : z  max 

z(i)  =  z_0  +  (i-l)*dz; 

end 

tic 

for  k=l:num_sys 

%sets  up  image  plane  grid 
dx  =  zeros  (z  max, 1) ; 
dy  =  zeros (z_max, 1) ; 
x  tot  =  zeros (x  max, z  max) ; 
y_tot  =  zeros  (y_max,  z_jnax)  ; 

for  q  =  1 : z_max 

%adjusts  image  plane  grid  magnification  based  on  distance  from 
%reference  plane  z_0 

dx(q)  =  dx_base  *  ( (q-1 ) *dz+z_0 )  /  z_0 ; 
dy(q)  =  dy_base  *  ( (q-1 ) *dz+z_0 )  /  z_0; 
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%Build  RMC  Filename  - 

% ' num_^sys-RMC_CL_coord-num_grids-grid_size-R  det-mask  file-mask  sep- 
N  samp-runtime-tau-epsRMC-x  max-y  max-dx  base-dy_base-x  expect-y  expect- 
z  max-z_0-dz .mat ' 

A  name  = 

[num2str (z (q) ) ,  , num2str (mask  sep (k)  )  ,  '  , num2str (min (t  mask (k) ) )  ,  '  , num2s 

tr(x  max),  ,num2str(y  max),  ' , num2str (dx (q) ) ,  ' , num2str (dy (q) ) , '  ,num2st 

r (x^expect) , '  ' , num2str (y_expect) ,  ; , num2str (grid  size(k)),'  ,num2str(N  sam 

p(k)),  , num2str (samp  diff  var(k) ) , '  ( 1 , num2str (RMC_CL  coord ( k, 1 )), '  ,num2s 

tr (RMC_CL_coord (k, 2) ) , ' ) ' , , num2str (mask_file (k, : ) ) ] ; 

%Read  RMC  filename  to  see  if  library  exists.  For  random  sampling 
%it  will  re-compute  A  matrix  for  new  pints  even  though  the  number 
%of  sampled  spoints  is  the  same. 

%checks  to  see  if  the  TF  is  sampled  randomly  or  at  a  constant 
%interval  and  recomputes  system  matrix  every  time  for  random 
if  var (dif f ( samp) ) <1 ; 

[A, x, dummy, dummy, y]  =  Read_RMCLib (A  name, Home^Address) ; 

else 

A=-l; 

end; 

if  (a  ==  -1) 

%Initialize  x  &  y 
x  =  zeros (x  max,l); 
y  =  zeros (y_max, 1) ; 

%****************************************  compute  RMC 
for  i=l:grid  size 

grid_dim(i)  =  (2*R_det  (k)  *  (i-1)  /  (grid__size  (k) -1)  ) 

end 

Q-'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 
o 

%***  Compute  grid  magnification  factor  'M' 

M  =  z  (q)  .  /  (z  (q)  -L  (k,  :  )  )  ; 

^★'^ryr'^ryr'^rVr'^ryr'^rVr'^ryr'^ryr'^rVr'^rVr'^ryr'^ryr'^rVr'^rVr'^rVr'^rVr'^ryr'^ryr 


grid  size 
-  R  det ( k) ; 


%Get  the  grid  patterns 
R  grids  =  R  det(k) . /M; 

R  slits (k)  =  R  det(k)  -  mask  min(k);  %  Sets  the  an  inner  radius 
as  a  contraint  on  the  slit  dimensions 
grids  = 

get_RMCgrids_v2 (a (k, : ) , w (k, : ) , c (k, : ) , M, grid_dim, R_slits (k) , slot_type (k, : ) , slo 
t  len ( :  ,  k)  )  ; 


O, 

o 

%  Calculate  the  system  matrix  'A' 

A  =  zeros (N  samp(k),x  max*y  max) ; 

h  =  waitbar ( 0 ,[' Pre-computing  A  -  ' , num2str ( ( k— 1 ) *z  max+q) , '  of 
,num2str(num  sys*z  max) ] ) ; 

for  i=l : x  max 
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x(i)  =  (i-round (x_max/2) - (x_expect/dx (q) ) )  *  dx (q) ; 

Sx  =  x(i)-RMC_CL  coord (k,l); 

for  j=l:y_max 

y(j)  =  (j -round (y_max/2 )- (y_expect/dy (q) ) )  *  dy(q); 

Sy  =  y ( j ) -RMC_CL_coord ( k, 2 ) ; 

A(:,  ( i-1 ) *y_max+ j )  = 

calcA  RMCgrids  v3 (grids, num  grids ( k) , grid  size (k) , Sx, Sy, z (q)  ,  L (k,  : )  ,  R  grids, g 
rid  dim, w ( k, : ) , M, N  samp ( k) , samp ( k, :), mask  sep ( k) , mask^coord) ; 

A(:, ( i-1 ) *y_max+ j )  =  (A ( : , (i-1) *y_max+j )  +  (1-A(:, (i- 

1) *y_max+j ) ) . *atten (k) ) . *RMC_Omega (R_det (k) , z (q) - 
max (L (k, : ) ) , sqrt (SxA2+SyA2) ) ./ (4*pi) ; 
end 

waitbar(i/x  max) 

end 

close (h) 

Write^RMCLib (A  name,A,x,x  max, dx (q) , y, y  max, dy (q) , Home_Address ) 

end 

if  (q  ==  1) 

A_^z  =  A; 
x_tot ( : , q)  =  x; 
y_tot ( : , q)  =  y; 

else 

A  z  =  [A  z ,  A]  ; 
x_tot ( : , q)  =  x; 
y_tot ( : , q)  =  y; 

end 

end 

A  =  A  z  ; 
clear  A  z; 

if  (k  ==1) 

A_tot  =  A; 

else 

A_tot  =  [A_tot;A] ; 

end 

clear  A 


end 

Write  RMCLib (filename, A  tot,x  tot,x  max, dx, y_tot, y  max, dy, Home  Address) 
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A.2  RMC  Toolbox.m 

%RMC  Toolbox 

%By:  Ben  Kowash 
%Date:  27  Jun  08 
%v :  1.0 

o, 

o 

%Description : 

%1 .  This  program  brings  together  many  of  the  RMC  analysis 
^concepts  I  have  developed  and  puts  them  into  a  single  package.  The 
%program  begins  by  asking  whether  the  input  data  will  be  provided 
%by  measured  data  or  by  simulated  data.  If  the  data  is  measured,  two  files 
%are  read  for  each  measurement  made:  (1)  Measured  data  for  each  measurement 
%and  (2)  System  configuration  for  each  measurement.  If  the  data  is  to  be 
%simulated,  a  pair  of  files  is  read  that  contains:  (1)  Location  and  activity, 
%of  all  sources  as  well  as  the  system  background.  (2)  The  system 
configuration 

%for  each  measurement.  To  avoid  confusion,  the  source  plane  is  used  to 
%define  the  coordinate  system  used  for  each  of  the  measurements 

%2 .  Once  the  problem  is  defined  a  system  model  is  generated  for  each  new 
%measurement .  The  user  defines  the  model  resolution  as  well  as  number  of 
%pixels  to  be  used, 
tic 

close  all 
clear  all 
clc 

num  meas  =  1;  %Number  of  independent  measurements 
simulate  =  'yes';  %Use  simulation  vs  measured  data 
filename  =  'Mask8-20cm' ;  %Name  of  case  to  run 
Home  Address  =  'F:\Lab  RMC! ' ; 

[A,x,x  max, dx, y, y_max, dy]  =  Read  RMCLib  (filename,  Home__Address)  ; 
plot  num  =  3; 

run  LL  =  'no';  %Source  estimation  with  log-likelihood  model 
run  fisher  =  'yes';  %Compute  Fisher  Information  map 
run  MLEM  =  'no';  %Image  reconstruction  using  ML-EM  algorithmE;\ 
iter  MLEM  =  100;  %#  of  MLEM  iterations 
iter_BS  =1;  %#  of  bootstrap  iterations 
MLEM  bkg  =  15;  %#  background  to  use  for  MLEM  algoritm 
MLEM  tol  =  le-5;  %Tolerance  for  MLEM  convergence 
z  max  =  1;  %#  of  z  planes  in  problem 


if  (strcmp (simulate, ' yes ' )  ==  1) 

sim  system  =  [Home_Address , ' RMCLib! Sim  System! ', filename, '  simsys . mat ' ] ; 
sim  source  =  [Home_Address , 'RMCLib\Sim  Src\ ', filename, '  simsrc .mat ' ] ; 
load ( sim_system) 
load (sim_source) 

else 

meas_system  = 

[Home_Address , ' RMCLib\Meas  System! ', filename, 1  meassys .mat ' ] ; 

meas  data  =  [Home  Address ,' RMCLib!Meas  data! ', filename,  measdata .mat ' ] ; 
load (meas_system) 
load (meas  data) 
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end 


Q-'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%*******************Set  Up  Data  Vector  T_sys************************* 

9r  ~k  9r  ~k  -k  -k  ~k  -k  ~k  9r  -k  9r  9r  9r  -k  9r  9r  9r  9r  9r  9r  9r  9r  9r  9r  9r  9r  -k  9r  9r  9r  9r  9r  9r  9r  9r  9r  9r  -k  9r  9r  9r  9r  -k  9r  9r  9r  9r  9r  9r  ~k  9r  9r  9r  9r  9r  9r  ~k  9r 

^* ************** QeD  simulated  system********************************* 

if  (strcmp (simulate, ' yes ' )  ==  1) 

T_sim  =  [ ] ; 
for  k=l:num_sys 

%***  Compute  RMC  grid  size 
for  i=l:grid  size(k) 

grid_dim(i)  =  (2*R_det  (k)  *  (i-1)  /  (grid__size  (k) -1)  )  -  R_det(k); 

end 

%***  Compute  grid  magnification  factor  'M' 

M  =  Sxyz (1,3)  ./  (Sxyz (1, 3) -L (k, : ) ) ; 

R  grids  =  R  det(k)  . /M; 

R  slit  =  R  det(k)  -  mask  min; 

%First  get  the  grid  patterns 
grids  = 

get_RMCgrids_v2 (a (k,  : ) , w (k,  : ) , c (k,  : ) , M, grid_dim, R__slit (k)  ,  slot_type (k,  : ) , slot 
_len ( : ,  k) ) ; % (a,  w, c, M, x, R_slit, slot_type, slot_len) 
grids ( 1 , : )  =  0 ; 

T_multsrc  =  []; 
for  i=l:num  src 
T_sngl  = 

calcA  RMCgrids  v3 (grids, num  grids ( k) , grid  size (k) , Sxyz (i, 1) +RMC_CL_coord (k, 1) 

, Sxyz ( i , 2 ) +RMC  CL_coord ( k, 2 ) , Sxyz ( i ,  3 )  ,  L ( k,  : )  ,  R  grids,  grid  dim, w ( k,  : ) , M, N  sam 
p ( k) , samp ( k, : ) , mask_sep ( k) , mask^coord) ; 

T_sngl  =  ( (T_sngl  +  (1-T_sngl) . *atten (k) ) . *  act(i)  .  * 
epsRMC (k) . *RMC_Omega (R_det (k) , Sxyz (i, 3) - 
max (L (k,  : ) )  ,  sqrt ( (Sxyz (i, 1) +RMC_CL_coord (k, 1) ) A2  + 

(Sxyz (i, 2) +RMC_CL_coord (k, 2) ) A2) ) . / (4*pi)  +  bkg ( i ) )  . *tau ( k) ; 

T  multsrc  =  [T  multsrc,T  sngl]; 

end 

T  sim  =  [T  sim;sum(T  multsrc, 2)]; 

T_sys  =  poissrnd (T_sim) ; 

end 

figure (plot_num) 
plot ( samp, T_sys , 'g*-') 
hold  on 

plot ( samp, T_sim,  ' b . -  '  ) 

%plot (T_sim. /4.3245e+003, 'b. ') 
plot  num  =  plot  num  +  1; 

clear  T  multsrc 
clear  T_sngl 

^*  **************  Load  Measured  Data*********************************** 
else 


T_sys  = 


zeros  (360, 1)  ; 
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o\°  o\° 


%T_sys  =  [ ]  ; 
num_sys=l  ; 
for  k=l:num_sys 

load (RMC_data_f iles ( k,  : ) )  ; 
%T_sys  =  [T_sys;  cnt ' ] ; 
T_sys  =  T_sys  +  cnt'; 

end 

num^sys  =  1; 

figure (plot  num) 

plot(T_sys, 'r.-') 

plot  num  =  plot_num  +  1; 


end 

ic*ic*icjcic*^*********'k-k'k'k'k-k'k'k'k-k'k'k'k-k'k-k'k-k'k'k'k-k'k'k'k-k'k-k'k-k'k-k'k-k'k'k'k-k'k-k'k'k-k'k'k'k-k'k-k'k 

-k-k'k-k-k-k-k-k-k-k'k-k-k-k'k-k-k-k'k'k-k-k-k-k-k-k'k-k-k-k-k-k-k-k'k-k-k-k-k-k-k-k'k-k-k'k-k-k-k-k-k-k-k'k-k-k-k-k-k'k-k'k-k-k-k-k-k'k 


%***************j_,oq-lJi]<ielihood 
if  (strcmp(run  LL, 'yes')  ==  1) 


Function****************************** 


act_ref  =  act;  %source  activity  [Bq] 
bkg_ref  =  bkg;  %background  [cps] 
tau_tot  =  tau(l,:); 
for  k=2:num_sys 

tau_tot  =  [tau_tot,  tau(k, :)]; 

end 


%Set  up  and  Run  the  bootstrap  algorithm 
if  (strcmp (simulate, ' yes ' )  ==  1) 

%Set  up  T_BS  for  simulated  data 
T_BS ( : , 1)  =  T_sys; 
for  j=2:iter_BS; 

T_BS ( : , j )  =  poissrnd (T_sim) ; 

end 

else 

%Set  up  T_BS  for  measured  data 

T  BS  =  Resample  RMC (RMC_raw_f lies ( k, : ) , RMC  home(k),iter  BS); 

end 


LL  =  zeros (x  max*y  max, iter  BS); 

h  =  waitbar(0, 'Mapping  RMC  Cost  Function'); 
for  i=l : iter  BS 
LL ( : , i)  = 

get_LLRMC(T  BS(:,i),A,x  max, y^max, epsRMC, act_ref , bkg_ref , tau  tot, num  sys,N  sa 
mp)  ; 


waitbar (i/iter_BS) 
end 

close (h) 
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o\°  o\° 


LL  mean  =  mean(LL,2); 

LL_im  =  reshape (LL  mean,x  max,y  max) ; 

figure (plot^num) 
subplot (4,3,1) 
imagesc (x,  y,  LL  im) 
axis  image 
axis  xy 
colorbar 

plot  num  =  plot  num+1; 

map  =  zeros (x  max*y_max,l) 

P_sum  =  sum (A, 1); 
for  i=l:x  max*y_max 
for  k=l : 360 

map(i)  =  map(i)  +  T_sys(k)  *  log (A (k, i) /P_sum (i) ) ; 

end 

end 

figure (plot  num) 

imagesc (x, y, reshape (map, x  max,y  max)) 

end 

%***************pj_s]ler  information**  *************************  ******** 
if  (strcmp(run  fisher, 'yes')  ==  1) 

act  ref  =  3.7E10  *  100E-6;  IReference  activity  used  for  Fisher 
Information  and  LL 

bkg^ref  =  8;  %Reference  background  used  for  Fisher  Information  and  LL 
[F_act, F_x, F_y, CRLB_act, CRLB_x, CRLB_y]  = 
get  FisherRMC (A, x, y, x  max,y  max,N  samp, act^ref , bkg  ref , tau, epsRMC, num  sys,R  d 
et, Sxyz (1,3) ) ; 

figure (plot^num) 
subplot (3,2,1) 

imagesc (x, y, reshape (F_act, x_max, y_max) ) 

colorbar 

axis  xy 

axis  image 

title ([' F_a_c_t  -  N  = 

' , num2str ( sum (N  samp) ) ] , 'fontsize  ,16, ' f ontweight ' , ' b ' ) 
subplot (3,2,3) 

imagesc (x, y, reshape (CRLB  act,x  max,y  max)) 

colorbar 

axis  xy 

axis  image 

title ([' var_a_c_t  -  N  = 

' , num2str ( sum (N  samp) ) ] , ' fontsize  ,16, ' f ontweight ' , ' b ' ) 
subplot (3,2,2) 

imagesc (x, y, reshape (F  x  +  F^y,x  max,y  max)) 

colorbar 

axis  xy 

axis  image 
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title ( [ ' F_A_B  -  N  = 

' , num2str ( sum (N  samp) )], 'fontsize',16, ' f ontweight '  , ' b ' ) 
subplot (3,2,4) 

imagesc (x, y, reshape (CRLB  x  +  CRLB  y,x  max, y_max) ) 

colorbar 

axis  xy 

axis  image 

title (['var  A  B  -  N  = 

' , num2str ( sum (N  samp) ) ] , 'fontsize  ,16, ' f ontweight '  , ' b ' ) 

%Plot  position  variance  with  large  value  at  zero  removed 
tmpA  =  CRLB  x; 
tmpB  =  CRLB  y; 

for  qq  =  1:5; 

loc(qq)  =  f ind (tmpA+tmpB  ==  max (tmpA+ tmpB) ) ; 
tmpA (loc (qq) )  =  0; 
tmpB (loc (qq) )  =  0; 
end; 

subplot (3,2, 6 ) 

imagesc (x, y, reshape (tmpA  +  tmpB,x  max, y^max) ) 

colorbar 

axis  xy 

axis  image 

%  title (['var  A  Bother  -  N  = 

'  ,  num2str  ( sum  (N__samp)  )  ]  ,  '  fontsize  ',16,  '  f ontweight '  ,  '  b  '  ) 
title ([' varAB  - 

N=  , num2str ( sum (N_samp) ) ]  ,  fontsize',16,  1 f ontweight ' ,  1 b '  )  ; 
subplot (3,2,5) 

tmpABsq  =  reshape (tmpA  +  tmpB,x_max,y  max) ; 
for  i=l : x  max 

for  j=l : y  max 

if  ( j==round (y_max/2) ) 

tmpAB(i)  =  tmpABsq ( j , i ) ; 

end 


end 


end 

semi logy (x, tmpAB, 'm. - ' ) 
plot  num  =  plot  num  +  1; 

stats  =  [max (CRLB  act);  min (CRLB  act);  mean (CRLB  act);  max (CRLB  x  + 

CRLB  y) ;  min (CRLB  x  +  CRLB  y) ;  mean (tmpA  +  tmpB)]; 

%save ( [ ' C : \Users\Ben\Documents\PhD 

%Research\RMC_Simulations\Fisher  Inf ormation\N= ' , num2str (N_samp) , '  - 

L= ' , num2str (mask^sep) ,  ' cm. mat ' ] ,  ' F  act '  ,  ' F  alpha ' ,  ' F  beta ' ,  ' var^act ' ,  ' var  alp 
ha ' ,  ' var^beta ' ,  ' x '  ,  ' y '  )  ; 

end 
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%****************ML-em  Image  Reconstruction  with 
Bootstrap************************** 
if  (strcmp(run  MLEM, 'yes')  ==  1) 

%  load ( ' C : \Axy_W.mat ' ) 

%  x  max  =  180; 

%  y^max  =  180; 

tau_tot  =  tau(l,:); 
for  k=2:num  sys 

tau_tot  =  [tau_tot,  tau(k, :)]; 

end 

%Set  up  and  Run  the  bootstrap  algorithm 
if  (strcmp (simulate, ' yes ' )  ==  1) 

%Set  up  T_BS  for  simulated  data 
T_BS ( : , 1)  =  T_sys; 
for  j=2:iter_BS; 

T_BS ( : , j )  =  poissrnd (T_sim) ; 

end 

else 

%Set  up  T_BS  for  measured  data 
T_BS  =  zeros (N_samp ( 1 ), iter_BS )  ; 
num_sys  =  1; 
for  k=l:num  sys 

T_BS  =  T_BS  + 

Resample  RMC (RMC_raw  f iles (k,  : ) , RMC  home (1), iter  BS); 

%T_BS  =  T_sys; 

end 


end 

T_size  =  size (T_sys, 1) ; 

%load ( ' C : \Extended-T .mat ' ) 

%b  =  T_bkg;  %ones (T  size,l)  .*  MLEM  bkg  .*  tau  tot '  ; 
b  =  ones (T  size,l)  .*  MLEM  bkg  .*  tau  tot'; 

lambda  =  ones (x  max*y_max*z  max, iter  BS);  %Initialize  the  guess  for 
lambda 

A  =  A. *epsRMC (1) *tau_tot (1) ; 

h  =  waitbar ( 0 ,' Performing  bootstrap  analysis'); 
for  i=l : iter  BS 

%Run  the  iterative  ML-EM 

a  =  A'  *  ones (Tsize, 1 ) ; 

lambda  old  =  lambda ( :,i); 
gate  =  1; 
q_MLEM  =  1; 

while  (gate==l)  | |  (q_MLEM  <  10) 
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lambda (:,i)  =  RMC_MLEM ( lambda (:, i )  ,  A,  T_BS(:,i),  b,  a) ; 
pct_chg_lam  =  (lambda (:,i)  -  lambda_old) . /lambda_old; 
lambda_old  =  lambda (:, i) ; 

if  ( sum ( (pct_chg_lam  >  MLEM  tol))  >  0  &&  (q_MLEM  <=  iter  MLEM) ) 
q_MLEM  =  q_MLEM+l; 

else 

gate  =  0; 
if  (q_MLEM  <  10) 
gate  =  1; 

q_MLEM  =  q_MLEM  +  1; 

end 

end 


end 


lambda (:,i)  =  lambda (:,i)  ./  (3.7E7);  %Converts  lambda  from  Bq  to  mCi 

k=l; 


%  [img  max(i,l),  img  max(i,2)]  =  max ( lambda (:, i )) ; 

%  img  max(i,3)  =  x (floor ((img  max(i,2)/x  max (k) ) ) ) ; 

%  if  (floor (img  max(i,2)/x  max (k) )  ==  img  max(i,2)/x  max (k) ) 

%  img_max(i,4)  =  y(y_max(k)); 

%  else 

%  img  max(i,4)  =  y(img  max(i,2)- 

(floor((img  max(i,2)/x  max(k))))*y  max(k)); 

%  end 

waitbar (i/iter_BS) 

end 

close (h) 

%  lambda  =  lambda  ./  (epsRMC (k) . *sum (sum (tau) ) ) ; 

lambda  mean  =  mean (lambda, 2) ; 
lambda_std  =  std (lambda, 1, 2) ; 
lambda_pkvar  =  lambda  mean  ./lambda  std; 

lambda  tmp  =  lambda  mean; 
for  i=l:x  max*y_max*z  max 
lambda_tmp ( i )  =  []; 
lam  mn(i)  =  mean (lambda  tmp) ; 
lambda_SNR ( i )  =  lambda  mean(i)  /  lam  mn(i); 
lambda  tmp  =  lambda  mean; 

end 


lambda_SNR  =  lambda^mean  ./  lambda_pkvar ; 


o, 

o 


BS_stats ( k, 1 ) 
BS_stats ( k,  2 ) 
BS_stats ( k, 3 ) 
BS  stats ( k, 4 ) 


mean (abs (img  max ( : , 3) ) ) ;  %x  mean 
std (abs (img_max ( : , 3) ) ) ;  %x  err 
mean ( abs ( img  max ( : , 4 ) ) ) ;  %y  mean 
std (abs (img_max ( : , 4) ) ) ;  %y  err 
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o\°  o\°  o\°  o\°  o\°  o\°  o\°  o\°  o\°  o\°  o\°  o\°  H-  o\°  o\°  C/3  o\°  o\°  f— 1  o\°  o\°  o\°  o\° 


%  BS_stats 

if  (k  ==  1) 
clear  k 

end 

for  q=l : z  max 

figure (plot_num) ; 

plot  num  =  plot  num  +  1; 

%  subplot  (2,2,1) 

lambda  im  =  reshape ( lambda  mean ( (q- 
l)*x  max*y__max+l :  q*x  max*y  max)  ,x  max,y  max); 
set  (gca,  ' YDir ' ,  ' normal ' ) ; 
imagesc (x ( : , q) , y ( : , q) , lambda_im) 
axis  xy 
axis  image 

title ( ' MLEM  Reconstruction'); 

colorbar 

hold  on 

contour (x ( : , q) , y ( : , q) , lambda_im) 

title ( ' MLEM  Reconstruction '  ,  'fontsize',16,  ' f ontweight '  ,  ' b ' ) 
xlabel ( ' X  [cm]  '  ,  'fontsize',16,  ' f ontweight ' ,  ' b ' ) 
ylabel('Y  [cm] ', 'fontsize',16, ' f ontweight ' , 'b ' ) 
subplot (2 , 2 , 3 : 4 ) ; 

hist (reshape (lambda  im, x  maxA2, 1) , 100) ; 


subplot (2,2,2) 

lambda  im  =  reshape ( lambda  mean((q- 
)  *x  max*y__max+l:q*x  max*y_max),x  max,y_max); 

sig  =  max (max ( lambda  im (1 : length (x) /2 , (length (x) /2 ) +1 : length (x) ))) ; 
noise  = 

sum ( sum ( lambda  im ( 1 : length (x) /2 ,( length (x) /2 ) +1 : length (x) )) )  - 

ig) / ( length (lambda  im ( 1 : length (x) / 2 , ( length (x) /2 ) +1 : length (x) ) ) A2-l ) ; 

SNR  =  sig  /  noise; 

mage sc (x ( ( length (x) /2 ) +1 : length (x) , q) , y ( 1 : length (y) /2 , q) , lambda_im ( ( length (x 
/2  )  +1 : length (x) , 1 : length (x) /2 )  ' ) 
axis  xy 
axis  image 

img  title  =  strcat('MLEM  Reconstructed  Image:  SNR= ' , num2str (SNR) ) ; 
title (img_title) ; 
colorbar 
hold  on 

contour (x ( : , q) , y ( : , q) , lambda_im ' ) 
title(img  title,  'fontsize',16,  ' f ontweight '  ,  'b') 
xlabel ( ' X  [cm]  ',  'fontsize',16,  ' f ontweight ' ,  ' b ' ) 
ylabel ( ' Y  [cm] ' , ' font size ',16, ' f ontweight ' , ' b ' ) 

subplot (2,2,2) 

lambda  imstd  =  reshape (lambda_pkvar ( (q- 
l)*x  max*y_max+l : q*x  max*y  max) ,x  max, y_max) ; 

%  imagesc (x(:,q) ,y(:,q) , lambda_imstd) 

%  axis  xy 

%  axis  image 

%  title (' Standard  Deviation  of  MLEM  Image') 

%  colorbar 
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%  subplot (2,2,3) 

%  lambda  imSNR  =  reshape ( lambda_SNR ( (q- 

l)*x  raax*y_raax+l : q*x  max*y_max),x  max, y_max) ; 

%  imagesc (x ( : , q) , y ( : , q) , lambda_imSNR) 

%  axis  xy 

%  axis  image 

%  title (' Signal  to  Noise  Ratio') 

%  colorbar 

%  if  (z  max  >  1) 

%  plot_num  =  plot  num  +  1; 

%  end 

O, 

o 

%  subplot (2,2,4) 

%  plot (x ( : , q) , lambda  im (21,  :)  ,  ' b . - ' ) 

end 

end 

toe 

A. 3  correct  for  num  sys.m 

%%This  eliminates  the  needs  to  replicate  the  same  values  repetitively  for 
%%each  system  parameter. 

function  [corrected_parameter ]  =  correct_for_num_sys (num_sys, sys_parameter) 
if  size (sys_parameter, 1)  ~=  num^sys 

for  i  =  l:num_^sys; 

corrected_parameter (i, : ) =sys_parameter ; 

end; 


else 


corrected_parameter  =  sys_parameter ; 


end 

A.4  Read  RMCLib.m 

function  [A  rd,x  rd,  xmax  rd,dx  rd, y_rd, ymax  rd,dy_rd]  = 

Read  RMCLib(lib  title, Home  Address) 

%By:  Ben  Kowash 
%Date:  17  Jun  08 

o, 

o 

%v  1.0 

%Description :  This  function  takes  in  data  about  a  system  run  and  determines 
%if  a  file  already  exists  in  the  library.  If  it  does  the  library  file  is 
%opened  and  returned.  Otherwise  an  empty  data  set  is  returned. 

lib_path  =  [Home  Address, ' RMCLibX System  Matrix  1\ ' ] ; 

lib_vol  =  [ lib_path, lib  title, '.mat']; 
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2) 


if  (exist(lib  vol, 'file')  == 
load ( lib_vol ) ; 

A_rd  =  A; 
x  rd  =  x; 
xmax  rd  =  x  max; 
dx_rd  =  dx ; 
y_rd  =  y; 
ymax^rd  =  y_max; 
dy_rd  =  dy; 

else 

A  rd  =  - 1 ; 
x_rd  =  - 1 ; 
xmax  rd  =  -1; 
dx_rd  =  - 1 ; 
y_rd  =  - 1 ; 
ymax  rd  =  -1; 
dy_rd  =  - 1 ; 


end 

A. 5  get_RMCgrids_v2.m 

function  [grids_out]  =  get_RMCgrids_v2 (a, w, c, M, x, R_slit, slot_type, slot_len) 

%  Title:  get  RMCgrids 
%  By:  Ben  Kowash 
%  Date:  20  May  08 
%  version :  2.0 

o, 

o 

%  Description:  This  function  is  used  to  set  up  the  various  masks  that  are 
%  used  to  model  a  rotational  modulated  collimation  detector.  The  final 
%  output  is  a  vector  called  'grids_out'  that  contains  the  raw  grid 
%  information  for  'N'  masks.  This  new  version  uses  heaviside  functions  to 
%  generate  the  masks. 

num  grids  =  size (M, 2); 
grid  size  =  size  (x, 2); 

num^slots  =  ceil (-c*2/a) ; %size (slot_len, 2) ; 
grids  out  =  zeros (grid  size, num  grids); 

for  i=l:num  grids 

for  k=l:num  slots 

if  strcmp (slot_type, ' trapazoid ' ) ; 

%comptues  the  secant  line  that  forms  the  slanted  ends  of  a 

trapaziod 

x_start  =  c (i) +a (i) * (k-1) ; 

x  end  =  c  (i) +w (i) +a (i) * (k-1) ; 

y_start  =  sqrt (R_slitA2-x_startA2 ) ; 

y  end  =  sqrt (R_slitA2-x_endA2 ) ; 

slope  =  (y_end-y_start) / (x_end-x_start) ; 

y_int  =  y_start  -  slope* (x_start) ; 
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trap_slot_len  =  slope* (x)  +  y_int; 

for  j=l:grid  size 

HI  =  heaviside (x (j )- (a (i) * (k-1) +w (i) +c (i) ) /M (i) ) ; 

H2  =  heaviside (x (j )- (a (i) * (k-1 ) +c (i) ) /M (i) )  ; 
grids_out ( j , i )  =  grids^out ( j , i )  +  ( H2 ' — 

HI ' ) *trap_slot_len ( j ) /M ( i ) ; 

end 

%  subplot (4,1,1) ; plot (x, gridsout ( : , 1 ) , ' b . - ' ) ; 

elseif  strcmp (slot_type, 'rounded'); 

round_slot_len=sqrt (R  slitA2-x. A2) ; 
for  j=l:grid  size 

HI  =  heaviside (x (j )- (a (i) * (k-1) +w (i) +c (i) ) /M (i) ) ; 

H2  =  heaviside (x (j )- (a (i) * (k-1 ) +c (i) ) /M (i) )  ; 
grids_out ( j , i )  =  grids_out ( j , i )  +  ( H2  '  — 

HI ' ) *round  slot  len(j)/M(i); 

end 

elseif  strcmp (slot_type, 'rectangle'); 

x_start  =  c  (i) +a (i) * (k-1) ; 
y_start  =  sqrt (R_slitA2-x_startA2 ) ; 

HI  =  heaviside  (x- (a  (i)  *  (k-1 ) +w  (i) +c  (i)  ) /M  (i)  )  ; 

H2  =  heaviside (x- (a ( i )*( k-1 ) +c ( i )) /M ( i )) ; 

grids_out ( : , i )  =  grids_out ( : , i )  +  (H2 ' -HI ' ) *y_start/M ( i ) ; 

elseif  strcmp (slot_type, 'specified'); 

HI  =  heaviside  (x- (a  (i)  *  (k-1 ) +w  (i) +c  (i)  ) /M  (i)  )  ; 

H2  =  heaviside(x-(a(i)*(k-l)+c(i))/M(i)); 

grids_out ( :  ,  i )  =  grids_out ( : , i )  +  (H2 ' -HI ' ) *0 . 5*slot_len (k) /M (i) 


end 


end 


end 

%  subplot (4,1,1) 

%  plot (x, grids_out ( : , 1 ) , ' b . - ' ) 
%  subplot (4,1,2) 

%  plot (x, grids_out ( : , 2 ) , ' b . - ' ) 
%  subplot (4,1,3) 

%  plot (x, gridsout ( : , 3 ) , ' b . - ' ) 
%  subplot (4,1,4) 

%  plot (x, grids  out ( : , 4 )  ,  ' b . - ' ) 
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A. 6  calc_RMCgrids_v3.m 

function  [A_out]  = 

calcA  RMCgrids  v3 (grids  in,num  grids, grid  size, Sx, Sy, Sz, L, R  grids , grid^dim, w, 
M, N_samp, samp, mask  sep,mask  coord) 

%  Title:  calcA  RMCgrids 
%  By:  Ben  Kowash 
%  Date:  12  Dec  07 

o, 

o 

%  Description:  This  function  computes  the  transmission  probability  that  is 
used 

%  by  the  ML-EM  algorithm  to  reconstruct  an  image  based  on  measured  RMC 
%  data.  Tprob  uses  the  project  and  shift  method  to  find  the  transmission 
%  factor  given  a  particular  source  location. 

%This  sets  the  algorithm  up  for  values  that  fall  outside  of  the  RMC  FOV 
T_prob  =  zeros (N^samp, 1) ; 
phi_crit  =  -1; 

FOV^RMC  =  0.5* (6 . 7/mask_sep) *Sz; 

FOVjmask  =  abs (w (2 ) / (M (2 ) * ( 1-M (2 ) ) ) ) ; 
if  (sqrt (SxA2+SyA2)  >  FOV_RMC) 

%Phi  crit  is  the  width  of  the  cone  of  acceptance  for  the  masks 
phi_crit  =  2*abs (90-acosd (FOVj_mask/sqrt (SxA2+SyA2) ) ) ; 
if  (Sy  ==  0) 

phi_phase  =  90-0 . 5*phi_crit; 

else 

phi_phase  =  atand (Sx/Sy) -0 . 5*phi_crit; 

end 


end 

%If  the  masks  are  symmetric,  k  can  go  from  1:180  and  the  second  half  is 
%just  the  same  pattern  repeated.  This  conditional  tests  the  masks  and  if 
%they  are  symmetric  speeds  up  the  calculation. 

if  ( sum (dif f (w) )  ==  0) 

num_pts  =  round (N_samp  /  2); 

else 

num_pts  =  N  samp; 

end 

%Compute  area  of  intersection  between  the  two  masks 

[grid_center , Sx_mov, Sy_mov]  =  rotate_source (num_grids , Sx, Sy, 1 , M, N_samp) ; 

%define  distance  between  detector  plane  mask  and  top  RMC  mask  in  x-y  plane 
d  =  sqrt ( (grid  center (num  grids , 1 ) -grid_center ( 1 , 1 )) A2  + 

(grid  center (num  grids , 2 ) -grid_center ( 1 , 2 ) ) A2 ) ; 

A  tot  =  (pi*R_yjrids (num  grids) A2);  %R  grids (num  grids) 

if  ( (R  grids (num  grids) +d)<=  R  grids (1)) 

A  int  =  A^tot; 

else 


%Define  angles  between  center  of  circles  and  points  of  intersection  on 
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%circles 

alpha  =  2*acos ( (R  grids (1)A2  +  dA2  -  R  grids (num  grids) A2)  / 

(2*R_grids (1) *d) ) ; 

beta  =  2*acos ( (R  grids (num  grids) A2  +  dA2  -  R  grids (1)A2)  / 

(2*R  grids (num  grids) *d) ) ; 

%Define  the  total  intersection  area 

A  int  =  0.5*alpha*R  grids(l)A2  -  0.5*R  grids ( 1 ) A2*sin (alpha)  + 

0.5*beta*R  grids (num  grids) A2  -  0.5*R  grids (num  grids) A2*sin (beta) ; 
end 

%This  loop  computes  the  transmission  probability  as  a  function  of  rotation 
%angle  k. 

for  k=l:num_pts 

[grid  center, Sx  mov,Sy_mov]  = 
rotate_source (num_grids , Sx, Sy, samp ( k) , M, N_samp) ; 
if  (phi_crit  >  0) 

if  ( (phi_phase  >=  0)  &&  (samp(k)  >  phi_phase  &&  samp(k)  < 

phi_phase+phi^crit)  | |  ((samp(k)  >  180+phi_phase  &&  samp(k)  < 
180+phi_phase+phi_crit) ) ) 

T_prob(k)  = 

shift  RMCgrids  v5 (grids  in ( : , 1 : 2 ) , 2 , grid_size, grid  center ( 1 : 2 ,:), R  grids(l:2) 
,grid  dim, Sx  mov, Sy  mov, Sz, samp (k) )  ./  A  tot; 

elseif  ( (phi_phase  <  0)  &&  (samp(k)  <  (phi_phase+phi_crit) )  | | 

(samp(k)  <  phi_phase+180  &&  samp(k)  >  phi_phase+180+phi  crit)  | |  ((samp(k)  > 
360+phi_phase  &&  samp(k)  <  360+phi_phase+phi_crit) ) ) 

T_prob ( k)  = 

shift  RMCgrids  v5 (grids  in ( : , 1 : 2 ) , 2 , grid_size, grid  center ( 1 : 2 ,:), R  grids(l:2) 
,grid  dim,  Sx  mov,  Sy  mov,  Sz,  samp  (k)  )  ./  A_^tot; 

end 


else 


%calculates 
T_prob(k)  = 

shift  RMCgrids  v5 (grids  in, num  grids, grid  size, grid_center, R  grids , grid_dim, S 
x_mov, Sy_mov, Sz, samp (k) )  ./  A_tot; 

T_prob(k)  =  T_prob(k)  + 

shift  RMCgrids  v5 (grids  in ( :  ,  1 : 2 )  ,  2  ,  grid  size, grid  center ( 1 : 2 ,:), R  grids(l,l: 
2),  grid  dim, Sx  mov, Sy_mov, Sz, k) * (A  tot-A  int) /A_tot A2 ; 

%  area(k)  =  mask  matrix  builder(R  grids , grid_center, mask  coord); 

end 


end 

%  save('F:\Lab  RMC\Results\Movies\mask  movie .mat Mask8  16  O', 'area'); 
%  T_prob  =  area./A_tot; 

%  plot (T_prob, ' r . - ' ) ; 
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%  pause; 


if  (num_pts  ==  round(N  samp/2)) 

T_prob (round (N  samp/2) +1:N  samp)  =  T_prob ( 1 : round (N  samp/2)); 

end 


A_out  =  T_prob; 

A. 7  Write  RMCLib.m 

function  Write  RMCLib(lib  title,A, x,x  max, dx, y, y_max, dy, Home  Address) 

%By:  Ben  Kowash 
%Date :  17  Jun  08 

o, 
o 

%v  1.0 

%Description :  This  function  takes  in  data  about  a  system  run  and  writes  a 
%RMC  system  library  based  on  the  input  information. 

lib_path  =  [Home^Address RMCLib\ System  Matrix/']; 

lib_vol  =  [ lib_path, lib  title, ' .mat']; 

save  ( lib_vol ,  '  A '  ,  '  x  '  ,  '  x  max  '  ,  '  dx  '  ,  '  y '  ,  '  y  max  '  ,  '  dy '  ) 

A. 7  shift_RMCgrids_v5.m 

function  [T_out]  = 

shift  RMCgrids  v5 (grids  in,num  grids, grid  size, grid  center, R  grids, grid  dim, x 

p,  yp, zp, pos _ plot ) 

%  Title:  shift  RMCgrids 
%  By:  Ben  Kowash 
%  Date:  12  Dec  07 
%  Updated:  2  May  08 
%  version:  4 

O, 

o 

%  Description:  This  function  takes  a  vector  of  grids  in  and  then  shifts  the 
%  grids  with  respect  to  the  first  grid  dependent  on  the  location  of  the 
%  source.  The  function  "get_pro j "  is  used  in  this  routine  to  calculate  the 
%  center  of  the  new  shifted  grid  so  that  the  new  position  can  be  computed. 

O, 

o 

%  In  this  updated  version  I  find  the  edges  of  the  slits  using  Matlab's  FIND 
%  command  and  then  rather  than  computer  the  solid  angle  I  just  find  the 
%  ratio  between  the  area  of  the  rectangles  and  the  entire  disk.  This  gives 
%  me  the  %  photons  transmitted. 

O, 

o 

%  Differences  from  v.3:  Instead  of  searching  the  entire  grid  of  Q1  = 

%  prod(masks),  I  now  circshift  Q1  and  take  the  difference.  This  then 
%  leaves  only  the  points  where  the  data  changes  from  zero  to  1 . 

O, 

o 

%  v.5  allows  the  function  to  compute  the  mask  area  of  trapazoid  slots.  This 
%  includes  finding  the  right  and  left  edge  of  every  slot  using  a  slightly 
modified  circshift 

%  method  from  previous  versions.  It  also  finds  the  corners  of  the  trapzoid 
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%  with  two  top  edge  (TE)  points  and  two  bottom  edge  points  (BE) . 

grid  overlap  =  ones (num  grids, 1); 
grid  =  zeros (grid^size, num  grids); 
grid(:,l)  =  grids^in ( : , 1 ) ; 

for  i=2:num  grids 

%Find  the  left  &  right  edges  to  ensure  the  grids  overlap 
LE  =  grid  center (i,l)  -  R  grids (i);  %Left  edge  of  grid  i 
RE  =  grid  center (i,l)  +  R  grids (i);  %Right  edge  of  grid  i 
if  (LE  >  R  grids (1))  | |  (RE  <  -R  grids (1)) 

grid_overlap ( i )  =  0; 

end 

if  (grid  overlap  ~=  0) 

%Shift  grid  ' i '  based  on  the  location  of  the  grid  ' i '  center  wrt  grid 

1 

shift  =  round (grid_center ( i , 1 )  /  (2*R  grids (1))  *  grid_size) ; 

if  (shift  >=  0) 

grid (shift+1 :grid_size, i)  =  grids  in(l:grid  size-shift, i) ; 
grid (1 : shift, i)  =  0; 
elseif  (shift  <  0) 

grid(l:grid  size-abs (shift) , i)  = 
grids  in (abs (shift) +1 :grid_size, i) ; 

grid (grid_size-abs (shift) :grid  size,i)  =  0; 

end 

end 


end 


Q1  =  prod (grid, 2 ) ; 

Q2  =  abs (Ql-circshift (Ql ,  1 ) )  ; 

%  figure  (3 )  ; 

%  subplot ( 1 , 2 , 1 ); plot (Q2 ,' b ', Ql ,' r '); hold  on; 

%  subplot ( 1 , 2 , 1 ); plot (Ql ,' r '); hold  off; 
first  =  find(Q2  ~=  0); 

Q2  =  round (Q2 . /min ( [Ql ( first ( 1 ) )  Q2 (first (length (first) ))])) ; 

mask_area  =  0; 

rect_edges  =  find(Q2  ~=  0); 

for  i  =  2 : 2 : length (rect_edges) ; 

rect_edges ( i )  =  rect_edges ( i )  -  1; 

end; 

rect  size  =  size  (rect  edges, 1) ; 

j  =  i; 

%  figure (4); 

p (1 : 5, 1 : 2 , 1 : rect_size-l)  =  0; 
while  j<=rect_size-l 

LE  =  grid_dim ( 1 , rect_edges ( j ) ) ; 

TE(1)  =  min (grid (rect_edges (j ),: )  +  grid_center ( : , 2 ) ' ) ; 
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BE(1)  =  max (grid_center ( : , 2 ) '  -  grid (rect_edges ( j ) , : ) ) ; 

%  if  (Q2 (rect_edges ( j+1) )  ~=  0) 

RE  =  grid_dim ( 1 , abs (rect_edges ( j +1 )  )  )  ; 

TE(2)  =  min (grid (rect_edges (j+1)  ,: )  +  gridcenter ( : , 2 ) ' ) ; 

BE (2)  =  max (grid_center ( : , 2 ) '  -  grid (rect_edges (j+1) ,:)) ; 

%  else 

%  LE  =  0; 

%  RE  =  0; 

%  TE  (  :  )  =  0 ; 

%  BE  (  : )  =  0 ; 

%  end 

p  (  : ,  : , j )  =  [LE  BE  (1)  ;  LE  TE  (1)  ;  RE  TE(2);  RE  BE(2);LE  BE ( 1 )  ] ; 

mask  area  =  mask  area  +  (RE-LE) * (min (TE) -max (BE) ) + . 5*abs (TE ( 1 ) 
TE (2) ) * (RE-LE) +. 5*abs (BE (1) -BE (2) ) * (RE-LE)  ; 

j  =  j+2; 


end 


%  figure (1) ;plot (Ql) ; legend (num2str (mask^area) ) ;hold  on; 

%  for  j  =  1 : rect_size-l ; 

%  subplot  (1, 2,2)  ; 

%  plot (p ( : , 1, j ) ,p ( : , 2, j ) , 'b' ) ;hold  on; 

%  end; 

%  hold  off; 

x=-5 : .01:5; 

%  if  num_^grids  >  2; 

%  plot (grid_dim, grid (:, 1) r '); axis ( [-5  5  -5  5]  )  ; 

%  hold  on 

%  plot (x, sqrt (Rgrids (1) A2- (x- 

grid  center ( 4  ,  1 )  )  . A2 ) +grid_center (4,2),  ' b ' ) ; 

%  plot (x, -sqrt (R_grids ( 1 ) A2- (x- 

gridcenter (4,1) ) . A2 ) +gr id_center (4,2), ' b ' ) ; 

%  plot (grid_dim, -grid ( :  ,  1 )  ,  ' r '  ) 

%  plot (grid_dim, grid ( : , 2 )  +  gridcenter (2 , 2 ) ) ; 

%  plot  (grid_dim,  grid__center  (2 , 2  )  -  grid(:,2)); 

%  plot (grid_dim, grid ( : , 3 )  +  grid  center ( 3 , 2 ) ,  ' g  '  )  ; 

%  plot (grid_dim, grid_center ( 3 , 2 )  -  grid ( : , 3 ) , ' g ' ) ; 

%  plot (grid_dim, grid ( : , 4 )  +  grid_center ( 4 , 2 ) , ' m ' ) ; 

%  plot (grid_dim, grid_center ( 4 , 2 )  -  grid ( : , 4 ) ,  ' m ' )  ; 

%  legend (num2str (mask  area)) 

%  pause (0.1) 

%  hold  off 

%  else 

%  figure  (2 ) ; 

%  plot (grid_dim, grid (:, 1)  ,' r '); axis ( [-5  5  -5  5]  )  ; 

%  hold  on 

%  plot (grid_dim, -grid ( : , 1 ) , ' r ' ) 

%  plot (grid_dim, grid ( : , 2 )  +  gridcenter (2 , 2 ) ) ; 
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%  plot (grid_dim, grid_center (2 , 2 )  -  grid(:,2)); 

%  %  plot (grid_dim, grid ( :  ,  3 )  +  grid_center ( 3 , 2 ) , ' g ' ) ; 

%  %  plot (grid_dim, grid_center ( 3 , 2  )  -  grid ( : , 3 ) ,  ' g ' )  ; 

%  %  plot (grid_dim, grid ( :  ,  4 )  +  grid_center ( 4 , 2 ) , ' y ' ) ; 

%  %  plot (grid_dim, grid_center ( 4 , 2  )  -  grid ( : , 4 ) ,  ' y  '  )  ; 

%  legend (num2str (mask  area)) 

%  pause  (0.1) 

%  hold  off 

%  end; 

if  (prod (grid^overlap)  ==  0) 

T_out  =  0 ; 

else 

T_out  =  mask_area; 

end 

A. 8  mask  matrix  builder.m 

function  [open  area]  =  mask  matrix  builder (R  grids, grid  center , mask^coord) 

%Mask  Matrix  Builder 
%18  Jan  2011 


r  =  R  grids ( 1 ) ; 
step  =.01; 

M  =  R  grids (1) ./R  grids; 

grid_center ( : , 1 )  =  grid_center(:,l)-grid_center(4,l); 
grid_center ( : , 2 )  =  grid_center ( : , 2 ) -grid_center ( 4 , 2 ) ; 

xl  =  -r:step:r; 

%  yl  =  sqrt (rA2-xl.A2) ; 

%  y2  =  -sqrt (r A2-xl . A2 ) ; 

O, 

o 

%  y3  =  -r:step:r; 

%  x2  =  sqrt (rA2-y3 . A2 )  ; 

%  x3  =  -sqrt (rA2-y3 . A2) ; 

theta_step  =  .01; 

theta  =  0:theta_step:2*pi+theta_step; 

[xy_coord ( : , 1 , 1 ) , xy_coord ( : , 2 , 1 ) ]  =  pol2cart (theta, r); 


%  xy^coord ( : , 1 , 1 )  =  cat (2,  cat (2, xl, xl) , cat (2, x2, x3) ) ; 

%  xy^coord  ( : ,  2 , 1 )  =  cat  (2,  cat  (2,  yl,  y2)  ,  cat  (2,  y3,  y3)  )  ; 

mask_coord  =  mask__coord  (1 :  4  ,  :  ,  :  )  ; 
slot_coord (:,:,:, 1 )  =  mask^coord; %*2 . 54 ; 

[~, ~, num^slots]  =  size (mask_coord) ; 

=  1:4; 

xy_coord ( : , 1 , i )  =  (xy_coord(:,l,l)./M(i))+grid_center(i,l); 
xy_coord ( : , 2 , i )  =  (xy_coord ( : , 2 , 1 ) . /M ( i ) ) +grid_center ( i , 2 ) ; 
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for  i 


slot_coord ( : , 1 , : , i )  =  ( slot_coord ( : , 1 , : , 1 ) . /M ( i ) ) +grid_center ( i , 1 ) ; 
slot_coord ( : , 2 , : , i )  =  ( slot_coord ( : , 2 , : , 1 ) . /M ( i ) ) +grid_center ( i , 2 ) ; 
xyloc ( : , 1 , i )  =  round (xy_coord (:, 1, i) /step)  +  length (xl ) +100 ; 
xyloc ( : , 2 , i )  =  length (xl ) +100  -  round(xy_coord(:,2,i)/step); 
slot_loc ( : , 1 , : , i)  =  round (slot_coord (:, 1, i) /step)  +  length (xl ) +100 
slot_loc ( : , 2 , : , i)  =  length (xl ) +100  -  round (slot_coord (:, 2, i) /step) 

end; 


%  figure (3) ;plot (xy_coord ( : , 1, 1) , xy^coord ( : , 2, 1) , ' r- ' ) ;hold  on; 
%  plot (xycoord ( : , 1 , 2 ) , xy_coord ( : , 2 , 2 ) , ' b— ' ) ; 

%  plot (xy  coord ( ; , 1 , 3 ) , xy_coord ( : , 2 , 3 )  ,  ' g- ' ) ; 

%  plot (xy_coord ( : , 1 , 4 ) , xy^coord ( :  ,  2 , 4 )  ,  ' y- ' ) ; 

gridl ( 1 : 2* length (xl ) +101 , 1 : 2* length (xl) +101) =0; 
grid2  =  gridl; 
grid3  =  grid2; 
grid4  =  grid3; 

%  plot (xyloc {: ,2 ,2) , xyloc ( : , 1 , 2 ) ) ; 
for  i  =  1: length (xy  loc) ; 

gridl (xy_loc ( i , 2 , 1 ) , xy_loc ( i , 1 , 1 ) ) =1 ; 
grid2 (xy_loc ( i , 2 , 2 ) , xy_loc ( i , 1 , 2 ) ) =1 ; 
grid3 (xy_loc ( i , 2 , 3 ) , xy_loc ( i ,  1 , 3 ) ) =1  ; 
grid4 (xy_loc ( i , 2 , 4 ) , xy_loc ( i , 1 , 4 ) ) =1 ; 

end; 


gridl 

=  imfill 

(gridl , 

' holes ' ) ; 

grid2 

=  imfill 

(grid2 , 

' holes ' ) ; 

grid3 

=  imfill 

(grid3 , 

' holes ' ) ; 

grid4 

=  imfill 

(grid4 , 

' holes ' ) ; 

det  = 

abs (1-gr 

id4 )  ; 

gridl  =  roipoly (gridl , xy_loc ( : , 1 , 1 ) , xy_loc ( : , 2 , 1 ) ) ; 
grid2  =  roipoly (grid2 , xy_loc ( : , 1 , 2 ) , xy_loc ( : , 2 , 2 ) ) ; 
grid3  =  roipoly (grid3 , xy_loc ( : , 1 , 3 ) , xy_loc ( : , 2 , 3 ) ) ; 
grid4  =  roipoly (grid4 , xy_loc ( : , 1 , 4 ) , xy_loc ( : , 2 , 4 ) ) ; 
det  =  abs ( l-grid4 ) ; 


slotsl  =  0; 
slots2  =  0; 
slots3  =  0; 
slots4  =  0; 


for  i  =  l:num  slots 
slotsl  =  slotsl 
slots2  =  slots2 
slots3  =  slots3 
slots4  =  slots4 

end; 


+  roipoly (gridl , slot^loc ( : , 1 , i , 1 ) , slot_loc ( : , 2 , i ,  1 ) 
+  roipoly (grid2 , slot_loc ( : , 1 , i , 2 ) , slot_loc ( :  ,  2  ,  i ,  2 ) 
+  roipoly (grid3 , slot_loc ( : , 1 , i , 3 ) , slot_loc ( :  ,  2  ,  i ,  3 ) 
+  roipoly (grid4 , slot_loc ( : , 1 , i , 4 ) , slot_loc ( :  ,  2  ,  i ,  4  ) 


gridl  =  gridl-slotsl ; 
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grid2  =  grid2-slots2 ; 
grid3  =  grid3-slots3; 
grid4  =  grid4-slots4 ; 

c  =  ceil ( (gridl+grid2+grid3+grid4 )  . /4  )  ; 
composite  =  c  +  det; 

composite  =  ceil (composite . /max (max (composite) ))  ; 
composite  =  abs (1-composite) ; 

open_area  =  length (nonzeros (composite) ) *stepA2 ; 
grey  =  - (c+det . * . 8) ; 

figure ( 4 ) ; image sc (grey (400:1500,200:1450) ) ; colormap ( ' gray ' ) ; axis 
off; title (open_area, ' Font Size ' , 24 ) ; 

O, 

o 

figure (5) ; imshow (composite (341:1186,341:1186),  ' InitialMagnif i cat ion '  ,  20) ;xlab 
el (open_area, ' Fonts ize ' , 24 ) ; 

%  pause  (.01)  ; 


A. 9  Resample  RMC.m 

function  [T  out]  =  Resample  RMC(raw  file, home, num^resample) 

%f  unction  [T_out]  =  Resample  RMC  (pos,  cnt,  home,  nurtures  ample) 

%Function  for  resampling  RMC  data  for  use  in  a  bootstrap  error  analysis 

O, 

o 

%  Ben  Kowash 
%  Date:  24  Apr  08 

O, 

o 

%  This  function  takes  the  raw  data  generated  by  the  Continuous  RMC 
%  Labview  routine,  and  resamples  it  to  generate  ' N'  independant 
%  realizations  of  the  same  modulation  profile.  The  output  is  a  length  PxN 
%  matrix,  where  P  is  the  #  of  independent  angles  of  rotation  for  the  RMC. 


%Open  the  raw  data  file 

load (raw_f ile) ; 

cnt  =  cnt ' ; 

pos  =  pos ' ; 

pos_start  =  home; 

%Reformat  the  position  data  array  into  angular  bins 
pos_size  =  size (pos, 1); 


o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 


%First  we  must  adddress  the  decoding  counter  overflow 
gate  =  1; 

stop_ovfl  =  pos__size; 
for  i=2:pos  size 

if  ((pos(i)  <  pos(i-l))  &&  (gate==l)) 

gate  =  0; 
start_ovfl  =  i; 
j  =start_ovf 1  +  1 ; 
while  ((pos(j)  >=  pos(j-l))) 

j  =  j+i; 

if  ( j==pos_size) 
break 

end 
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o\°  o\° 


overflow 


end 

%The  following  conditional  is  used  because  sometimes  the 


O, 

o 

and 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 


%counter  decrements  on  the  final  value...  i.e.  pos(i)  =  6300 

%pos(i-l)  =  6400.  Without  this  conditional  the  overflow  counter 
%would  be  off  by  one  position, 
if  ((pos(j)  >  3000)  &&  j<pos_size) 

j  =  j+i; 

end 

q, 

o 

stop_ovfl  =  j-1; 

pos ( start  ovf 1 : stop  ovf 1 )  =  65535  +  pos (start_ovf 1 : stop_ovf 1)  ; 


end 

if  ( i>stop_ovf 1 ) 
gate  =  1; 

end 


end 


for  i=l:pos_size 

if  (pos(i)  >  pos_start) 

pos(i)  =  pos(i)  -  pos_start; 

else 

pos(i)  =  72000  -  pos_start  +  pos(i); 

end 


end 

ang_pos  =  360  *  pos  /  (18000  *  4); 


%Count  #  of  total  rotations  in  file 
num^rot  =  1; 
for  i=2:pos_size 


if  (ang_pos(i)  <  ang_pos ( i — 1 ) ) 
num  rot  =  num^rot  +  1; 

end 

end 


T_rot  =  zeros (360, num^rot) ; 
gap  =  zeros (pos_size-l , 1 ) ; 

%Put  the  data  into  Q  independent  sampling  vectors.  Each  vector  here 

%represents  one  complete  rotation  of  the  RMC . 

k=l; 

for  i=2:pos_size 

incr  =  round (ang_pos (i) ) ; 
if  (incr==0) 

incr  =  360; 


end 
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prev  =  round (ang_pos (i-1) ) ; 
if  (prev==0) 

prev  =  360; 

end 

gap (i-1)  =  incr-prev; 
if  (gap(i-l)  <=2) 

T_rot ( incr , k)  =  T_rot ( incr , k)  +  (cnt (i) -cnt  (i-1) )  ; 

else 

for  j =prev : prev+gap (i-1) 

T_rot(j,k)  =  T_rot(j,k)  +  (cnt (i) -cnt  (i-1) ) /gap (i-1) ; 

end 

end 

if  (ang_pos(i)  <  ang_pos (i-1) ) 
k  =  k  +  1  ; 

end 


end 

%Error  check  to  make  sure  all  vectors  are  statistically  similar  --  this 
%is  important  if  the  RMC  only  completed  a  partial  rotation  on  the  final 
%pass 

T_sum  =  sum(T^rot) ; 

T  mean  =  mean (T  sum) ; 

T_std  =  std(T_sum); 

gate  =  0; 

while  (gate  ==  0) 

T  min  =  min (T_sum) ; 

min_pos  =  find(T_sum  ==  T  min); 

if  (T  min  <  T  mean  -  3*T  std) 

T_rot ( : , min_pos )  =  []; 

T_sum (min_pos )  =  []; 

T  mean  =  mean (T  sum) ; 

T_std  =  std(T_sum); 

else 

gate  =  1; 

end 


end 

num^rot  =  size (T_rot, 2 ) ; 

T_out  =  zeros(360,num_resample); 
T_out(:,l)  =  sum (T_rot, 2 ) ; 

for  i=2 : num_resample 

f  =  ceil (num  rot  . *  rand(num  rot,l)); 
T  out(:,i)  =  sum(T  rot(:,f),2); 


end 


end 
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A.  10  get  LLRMC.m 

function  [LL]  = 

get_LLRMC (T_in, A  in,x  max,y  max, epsRMC, act , bkg, tau, num^sys , N  samp) 

%Compute  the  log-liklihood  ratio  for  the  source  and  for  no  source. 
yn_bar  =  zeros (size (T_in, 2) , 1) ; 

tic 

for  i=l:x  max*y_max 

yn_bar  =  (A_in ( : , i ) . *epsRMC ( 1 ) . *act  +  bkg).*tau'; 

LL(i)  =  sum(T  in  .  *  log(yn  bar)  -  yn  bar); 


end 

A.  11  getFisherRMC.m 

function  [F  act,F  x,F_y,CRLB  act,CRLB  x,CRLB_y]  = 

get_FisherRMC  (A  in,x  in,y_in,x  max,y  max,N  samp, act, bkg, tau, epsRMC, num  sys,R 
det, Sz) 

N_samp_holder  =  1;  %A  Place  holder  to  keep  track  of  position  in  matrix  A 
tau_tot  =  tau(l,:); 
for  k=l:num__sys 

N  samp  holder  =  N_samp  holder  +  N  samp(k); 
if  (k>l ) 

tau_tot  =  [tau_tot, tau (k, : ) ] ; 

end 


end 

N  samp  tot  =sum (N_samp) ; 

d  alpha  =  zeros (N^samp  tot,x  max*y_max) ; 
d  beta  =  zeros (N^samp  tot,x  max*y__max)  ; 

N_samp_tot  =  360; 
for  k=l:N  samp  tot 

A  sub  =  reshape (A^in ( k, :), x  max,y  max) ; 

%f igure ( 1 ) 

%subplot (3,1,1) 

%imagesc (A_sub) 

%colorbar 
%pause (0.01) 

[d_aSub, d_bSub]  =  gradient (A_sub) ; 
d  alpha (k,  :)  =  reshape (d  aSub,l,x  max*y_max) ; 
d_beta(k, :)  =  reshape (d_bSub, 1 , x_max*y_max) ; 
%subplot (3,1,2) 

%imagesc (d_aSub) 

%colorbar 
%subplot (3,1,3) 

%imagesc (d_bSub) 

%colorbar 
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end 


clear  A_sub 

h  =  waitbar ( 0 Computing  Fisher  Information'); 
for  i=l : (x  max*y_max) 

pn_vec  =  zeros (1,3); 

F  =  zeros (3,3) ; 

Q  =  zeros (3,3) ; 

for  k=l:N  samp  tot 

if  (A  in ( k, i )  <  0 ) 

A_in ( k, i )  =  0  ; 

end 

pn_vec(l,l)  =  (tau(k)  *  epsRMC(l)  *  A_in(k,i)); 
pn_vec(l,2)  =  (tau(k)  *  act  *  epsRMC(l)  *  d_alpha ( k, i ) ) ; 
pn_vec(l,3)  =  (tau(k)  *  act  *  epsRMC(l)  *  d_beta (k, i) ) ; 

Q(l,l)  =  tau(k)  *  (act  *  epsRMC(l)  *  A_in(k,i)  +  bkg)  ; 

Q (2 , 2 )  =  Q (1, 1)  ; 

Q ( 3 , 3 )  =  Q (1, 1)  ; 

F  =  F  +  pn^vec '  *  pn  vec  /  Q; 
if  sum(F)  <  0 
pause 

end 

end 

F_act(i)  =  F ( 1 , 1 ) ; 

F_x(i)  =  F  (2 , 2  )  ; 

F_y(i)  =  F ( 3 , 3 ) ; 
sys_var  =  FA-1; 

CRLB_act(i)  =  sys_var ( 1 , 1 ) ; 

CRLB  x(i)  =  sys  var(2,2); 

CRLB_y(i)  =  sys  var ( 3 , 3 )  ; 

waitbar (i/ (x  max*y  max) ) 

end 

close (h) 

A.  12  RMC  MLEM.m 

function  lam  new  =  RMC_MLEM ( lam_old.  A,  y,  b,  a) 

%  Title:  RMC  ML-EM  routine 

O, 

o 

%  By:  Ben  Kowash 
%  Date:  15  Oct  07 

o, 

o 

%  Description:  This  routine  takes  input  data  and  computes  the  estimate 
%  lambda  new  based  on  the  properties  of  the  system.  The  method  used  is  a 
%  maximum-liklihood  expectation-maximum  iterative  algorithm. 
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lam  new 


lam  old  .  *  (A'  * 


(y  ./  (A 


lam  old  +  b)  ) )  .  /  a; 


A.  13  RMC  Omega.m 

function  g  =  RMC_Omega (R, z , rho) 

%  Ben  Kowash 
%  24  Sep  07 

%Description :  This  function  returns  the  solid  angle  of  a  point  source 
subtended  by  a 

%circular  disk  (i.e.  RMC  front  mask/detector) .  The  method  used  in  this 
%algorithm  is  provided  by  S.  Tryka  in  the  paper  "A  Method  for  Calculating 
%the  Average  Solid  Angle  Subtended  by  a  Circular  Disk  from  Uniformly 
%Distributed  Points  Within  a  Coaxial  Circular  Plane",  Review  of  Scientific 
%Instruments ,  Vol  70,  #10,  pg  3915-3920,  Oct  1999. 
switch  rho>=0 

case  rho==0 

g  =  2 . *pi . * (1-z . /sqrt (R. A2  +  z . A2) )  ; 
case  rho<R 

%Calc  parameters  for  elliptic  integrals  of  1st  &  3rd  kind 
Pp  sq  =  -4.*R.*rho  ./  ( (R-rho) . A2  +  z.A2); 
qp  =  -4.*R.*rho  /  (R-rho) .  A2; 

%Calc  elliptic  integral  of  1st  kind 

fun_K  =  @ (delta) (1./ (1-Pp_sq. *sin  ( delta) .A2) ,A0.5) ; 

K_Pp  =  quad ( fun_K, 0 , pi/2 , IE- 6) ; 

%Calc  elliptic  integral  of  3rd  kind 
fun_PI  =  @ (delta) (1./ ( (1-qp. *sin ( delta) .A2) . * (1- 
Pp_sq. *sin( delta)  .A2)  ,A0.5) ) ; 

PI  =  quad ( fun_PI , 0 , pi/2  ,  IE-6 )  ; 

%Calc  Omega 

g  =  (2.*pi)  -  (2 . *z . /sqrt ( (R-rho)  . A2  +  z . A2)  .  *  (K  Pp+ (R+rho)  . / (R- 
rho)  . *PI ) )  ; 

case  rho==R 


Ppr__sq  =  -4  .  *R .  A2  ./  z.A2; 

%Calc  elliptic  integral  of  1st  kind 

fun_K  =  @ (delta) (1./ (1-Ppr_sq. *sin( delta) .A2) ,A0.5) ; 

K_Ppr  =  quad ( fun_K, 0 , pi/2 , IE-6 ) ; 

%Calc  Omega 
g  =  (pi  -  2.*K  Ppr) ; 

case  rho>R 

%Calc  parameters  for  elliptic  integrals  of  1st  &  3rd  kind 
Pp_sq  =  -4.*R.*rho  /  ((R-rho).A2  +  z.A2); 
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qp  =  -4.*R.*rho  ./  (R-rho)  .  A2; 

%Calc  elliptic  integral  of  1st  kind 

fun_K  =  @ (delta) (1./ (1-Pp_sq. *sin  ( delta) .  A2) .  A0 . 5) ; 
K_Pp  =  quad ( fun_K, 0 , pi/2 , IE- 6) ; 

%Calc  elliptic  integral  of  3rd  kind 

fun  =  @  (delta)  (1 .  /  (  (1-qp.  *sin(  delta).  A2)  . *  (1~ 

Pp_sq. *sin (delta) .A2) . A0 . 5) ) ; 

PI  =  quadl (fun, 0, pi/2,  IE-6)  ; 

%Calc  Omega 
gl  =  -  2  .  *  z ; 

g2  =  sqrt ( (R-rho) . A2+z . A2 ) ; 

g3  =  (K  Pp+PI . * (R+rho) . / (R-rho) ) ; 

g  =  gl.*g3./g2; 

otherwise 

g  =  -100; 


end 

end 

A.  14  RMC_raw.m 

function  [T  dat]  =  RMC_raw (raw_file, RMC_file, home, plot_val) 

%Continuous  RMC  Raw  Data  Converter 

O, 

o 

%  Ben  Kowash 
%  Date:  8  Aug  07 

o, 

o 

%  This  subroutine  takes  the  raw  data  generated  by  the  Continuous  RMC 
%  Labview  routine,  and  converts  it  into  properly  binned  RMC  data  that  can 
%  be  processed  in  an  appropriate  routine  (ie  MLEM,  detector  algorithm,  etc) 

%Open  the  raw  data  file 

load (raw_f ile) ; 

cnt  =  cnt ' ; 

pos  =  pos ' ; 

pos_start  =  home; 

%Reformat  the  position  data  array  into  angular  bins 
pos_size  =  size  (pos, 1); 


for  i=l:pos_size 


if  (pos(i) 
pos (i) 

else 

pos (i) 

end 


>  pos_start) 

=  pos(i)  -  pos_start; 

=  72000  -  pos_start  +  pos(i); 


114 


end 


ang_pos  =  360  *  pos  /  (18000  *  4); 

if  (plotjval  ==  ' T '  ) 
figure  ( 1 ) 

plot (ang_pos ( : ) , ' g . ' ) 

end 

T_dat  =  zeros  (360, 1)  ; 
gap  =  zeros (pos_size-l, 1)  ; 
for  i=2:pos_size 

incr  =  round (angjpos (i) ) ; 
if  (incr==0) 

incr  =  360; 

end 

prev  =  round (ang_pos ( i — 1 ) ) ; 
if  (prev==0) 

prev  =  360; 

end 

gap(i-l)  =  incr-prev; 
if  (gap ( i — 1 )  <=2) 

T_dat ( incr ,  1 )  =  Tjdat ( incr , 1 )  +  (cnt (i) -cnt ( i — 1 ) ) ; 

else 

for  j=prev:prev+gap ( i  —  1 ) 

T_dat(j,l)  =  T_dat(j,l)  +  (cnt (i) -cnt ( i — 1 ) ) /gap ( i — 1 ) 

end 

end 


end 

cnt  =  Tjdat ' ; 

if  (exist ( [RMC_file] ) ==0) 

save ( [RMC_f ile] , ' cnt ; , ' pos ' ) 

['New  File  Written  =  , RMCJfile] 

else 

[  ' The  file  " ' , RMC  file,  ' "  already  exists !  ' ] 

end 

if  (plotjval  ==  ' T ' ) 
figure  (2 ) 
plot (T_dat,  ' b-o  '  ) 

end 

A.  15  rotate  source.m 

function  [grid  center,  Sx  mov,  Syjov]  = 
rotate_source (numjgrids , Sx, Sy, angular_pos , M, N_samp) 

%  Title:  rotate_source 
%  By:  Ben  Kowash 

115 


Date:  12  Dec  07 


o, 

o 

%  Description:  This  function  takes  in  the  position  of  the  RMC  system  and 
%  the  location  of  the  source.  The  location  of  the  projected  grid  center 
%  for  each  grid  is  then  output. 


%Grid  center  contains  the  x  &  y  coordinates  of  the  projected  mask  center 
grid  center  =  zeros (num  grids, 2); 

omega  =  2 *pi*angular_pos  /  N  samp; 


%Def  j 

if  ( 


:ine  the 
(Sx  >  0) 
Sx  mov  = 
Sy_mov  = 
elseif  ( (Sx 
Sx  mov  = 
Sy  mov  = 
elseif  ( (Sx 
Sx  mov  = 
Sy  mov  = 
elseif  ( (Sx 
Sx  mov  = 
Sy  mov  = 
elseif  (Sx  = 
Sx  mov  = 
Sy  mov  = 
elseif  (Sx  = 
Sx  mov  = 
Sy  mov  = 

end 


source  position 
&&  (Sy  >=  0) )  % 
^  sqrt (SxA2+SyA2) 

^  sqrt (SxA2+SyA2) 

<  0)  &&  (Sy  >=  0 
^  sqrt (SxA2+SyA2) 

^  sqrt (SxA2+SyA2) 

<  0)  &&  (Sy  <  0) 

:  sqrt ( SxA2  +  SyA2  ) 

:  sqrt ( SxA2  +  SyA2  ) 
>  0)  &&  (Sy  <  0) 

:  sqrt ( SxA2+SyA2 ) 

^  sqrt (SxA2+SyA2) 
;=  0)  &&  (Sy  >=0) 

^  sqrt (SxA2+SyA2) 

^  sqrt (SxA2+SyA2) 
:=  0)  &&  (Sy  <  0) 

^  sqrt (SxA2+SyA2) 

^  sqrt (SxA2+SyA2) 


for 

For 

* 


the  k-th  time  step 
src  in  1st  quadrant 
(cos (omega+atan (Sy/Sx) ) ) ; 
(sin (omega+atan (Sy/Sx) ) )  ; 
%For  src  in  2nd  quadrant 
(cos (omega+atan (Sy/Sx)  + 
(sin (omega+atan (Sy/Sx)  + 
For  src  in  3rd  quadrant 
(cos (omega+atan (Sy/Sx)  + 
(sin (omega+atan (Sy/Sx)  + 
For  src  in  4th  quadrant 
(cos (omega+atan (Sy/Sx) ) ) ; 
(sin (omega+atan (Sy/Sx) ) ) ; 


*  (cos (omega+pi/2 ) ) ; 

*  (sin (omega+pi/2) ) ; 

*  (cos (omega-pi/2 )) ; 

*  (sin (omega-pi/2)); 


PD 

Pi) 

Pi) 

Pi) 


%Next  find  new  location  of  A2center 
for  i=2:num  grids 


[grid_center ( i , 1 ) , grid_center ( i , 2 ) ]  = 
backpro j_RMC ( Sx_mov, Sy  mov, 0,0, M ( i ) ,0,0,0) ; 


end 


A.  16  ssim.m 

%Image  Quality  Assessment 

%This  will  compute  a  Structual  Similarity  Based  Image  Quality  Assessment 
%Index  (SSIM)  for  an  image  compared  to  a  truth  image. 

clear  all; 
close  all; 
clc; 

%Build  truth  matrix 
bk  =  15;  %cps 
src  =  .330;  %mCi 
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x  loc  =  10;  %x  pixel  location 
y_loc  =  -16;  %y  pixel  location 

num_pix  =  40;  %number  of  pixels  in  row  or  column 
%  loc  =  ( (sqrt (num_pix) /2) + (x_loc) - 

1) *sqrt (num_pix) + ( (sqrt (nurnjpix) /2) + (y_loc) ) ;  %find  location  in  array 
loc  =  [x_loc  num_pix+y^loc] ; 

truth ( 1 : num_pix, 1 : num_pix) =bk/ ( 3 . 7* 10 A7  )  ;  %sets  array  to  bk  in  mCi 
truth (loc (1) , loc (2) ) =src; 

for  j  =  1:8; 

path=' J:\Lab  RMC\Results\Ex  lOOmlem  100bs\'; 
file  =  [ 'Mask' , num2str ( j ) ,  -20cm' ] ; 

for  i  =  1:11; 

%load  image  data  from  figure  for  comparison 

if  i  ==  11; 

time  =  600; 

else 

time  =  30*i; 

end; 

name= [path, file, num2str (time) , ' s . fig ' ] ; 
load (name, ' -mat ' ) ; 
img  = 

reshape (hgS_070000 . children (1,1) . children (1,1) . properties . CData, 1 , num_pix) ; %c 
onverts  fig  data  into  array 

%Sets  parameters  for  SSIM 
mu  truth  =  mean (truth); 
mu  img  =  mean (img); 
var_truth  =  var (truth); 
var_img  =  var (img); 
sig_truth  =  sqrt (var_truth) ; 
sig_img  =  sqrt (var  img)  ; 
cov__mat  =  cov  (truth,  img)  ; 
cov_truth_img  =  cov  mat (1,2); 

alpha  =  1;  %These  are  using  to  adjust  importances  between  luminance, 
constrast,  and  structure 
beta  =  1; 
gamma  =  1 ; 

cl  =  0;  %These  constants  compensate  for  instabilities  for  near 

identical  differnces  in  the  denominator 
c2  =  0; 
c3  =  0; 

1  =  (2*mu  truth*mu  img+cl)/ (mu  truthA2+mu  imgA2+cl);  %luminance 
c  =  (2*sig  truth*sig  img+c2)/ (var  truth+var  img+c2);  %contrast 
s  =  (cov_truth_img+c3) / (sig_truth*sig_img+c3) ;  %structure 

SSIM(i,j)  =  lAalpha*cAbeta*s Agamma; 
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end; 


end; 

semi logy (cat (2 , 30:30:300, 600), SSIM, ' legend ( ' Mask  1 ' , ' Mask  2 ' , ' Mask 
3 ' , ' Mask  4 ' , ' Mask  5 ' , ' Mask  6 ' , ' Mask  7 ' , ' Mask  8 ' ) ; 

A.  17  backproj  RMC.m 

function  [x^center , y_center , w_pro j , a_pro j , c_pro j ]  = 
backproj  RMC(Sx  mov, Sy  mov, Ax_pr, Ay_pr, M, w, a, c) 

%  Function  backproj  RMC 

O, 

o 

%  By:  Ben  Kowash 
%  Date:  24  Oct  07 
%  Technique  by:  Ziya  Akcasu 
%  Version :  1.0 

O, 

o 

%  Description:  This  function  first  takes  a  mask  located  at  z=0  (i.e.  a 
%  mask  located  at  the  detector  plane)  and  backprojects  it  up  to  the  plane 
%  of  the  first  mask.  These  coordinates  are  then  used  to  compare  against 
%  the  first  mask  to  determine  the  total  detector  area  uncovered  as  the 
%  mask  system  rotates. 

%Compute  the  location  of  the  projected  mask  center 
%x_center  =  Sx  mov  +  (Sx  mov-Ax_pr) * (Az-Sz) / (Sz); 
x_center  =  Sx  mov  -  (Sx  mov-Ax_pr) /M; 

%y_center  =  Sy  mov  +  (Sy_mov-Ay_pr) * (Az-Sz) / (Sz); 
y_center  =  Sy  mov  -  (Sy  mov-Ay_pr) /M; 

%Compute  the  magnitude  of  the  projected  w. 
w_proj  =  (w/M) ; 

%Compute  the  magnitude  of  the  projected  a. 
a_pro j  =  (a/M) ; 

%Computer  the  magnitude  of  the  projected  c. 
c_proj  =  (c/M) ; 


A.  18  isEven.m 

function  [TF]  =  isEven(val) 

%By:  Ben  Kowash 
%Date :  26  Feb  08 

o, 

o 

%Description :  This  function  takes  in  a  value  and  determines  whether  that 
%value  is  odd  or  even.  If  the  value  is  even  a  1  is  returned.  Otherwise  a  0 
%is  returned  for  false. 

for  i=l : size (val) 

if  ( (mod (val ( i ) , 2 ) )  ==  0) 

TF ( i )  =  1; 
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else 


0; 


TF ( i )  = 

end 

end 


A.  19  Splice.m 

function  []  =  Splice (RMC_file, home, desired  time) 

%This  allows  for  a  long  RMC  measurement  to  be  spliced  into  shorter 
%measurements  to  look  at  the  time  dependency  of  reconstruction.  It 
%requires  that  the  RMC  file  index  be  a  three  digit  collected  time 
%identifier.  eg.  'C:\RMCout\RMC_450  -  Mask5-20cm.mat '  would  be  for  a  450 
%second  measurement  with  Mask5-20cm  being  the  title  entered  into  the  RMC 
%Continuous  VI  in  labview  or  from  the  RMC  command  txt  file.  It  also  usese 
%the  raw  data  file  associated  with  the  collection  time  index. 


collected  time  =  str2double (RMC_file (23 : 25) ) ; 

RMC_file  name  =  RMC  file (26 : length (RMC  file) ) ; 

raw_file  =  ['F:\Lab  RMC \RMCout\ data  num2str (collected  time), ' .mat']; 
load (raw_f ile) ; 

%Pulls  the  data  taken  during  the  desired  time 

cnt  =  cnt (1 : round (length (cnt) * (desired_time/collected_time) ))  ; 
pos  =  pos (1 : round (length (pos) * (desired_time/collected_time) )) ; 

new^raw  file  =  ['F:\Lab  RMC\RMCout\data  num2str (desired  time) ,' .mat '] ; 
new  RMC_file  =  ['F:\Lab 

RMC\RMCout\RMC  num2str (desired  time), RMC  filename]; 
save ( [new_raw_f ile] , ' cnt ' , ' pos ' ) ; 

RMC_raw (new_raw_f ile, new^RMC  file, home, 'F') 

A.20  Time  Detect.m 

function  [t  mean, t  std]  =  Time  Detect (cnt  bk, cnt  s  bk, collection  time) 

%This  computes  time  to  detect  based  on  one  activity  measurement  for  a 
%specific  channel  window  in  the  MCA  base  on  the  resolution  of  the  detector 
%for  a  specific  energy  level.  A  detection  threshold  of  photo-peak  height 
%above  background  is  used. 

%  clear  all; 

%  c  1  c  ; 

%  close  all; 

%standard  deviation  for  counts 
cnt_s_bk  =  sum(diff  (cnt_s_Jok)  )  ; 
cnt_s_bk_std  =  sqrt (cnt_s_bk) ; 
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cnt_bk_std  =  sqrt (cnt_bk) ; 


%activi ty 

s^bk  mean  =  cnt  s  bk/collection  time; 
s_bk_std  =  cnt_s_bk_std/collection_time; 
bk  mean  =  cnt  bk/collection  time; 
bk  std  =  cnt  bk  std/collection  time; 

src_mean  =  s_bk  mean  -  bk  mean; 
src_std  =  sqrt ( s_bk_stdA2  +  bk_stdA2); 

det  res  =7.5;  %  detector  resolution  in  percent  at  specific  energy  level, 
122.10  keV  for  Nal 
energy  =  122.10;  %  [keV] 

FWHM  =  (det  res/100) *energy; 
peak_std  =  FWHM/ (2*sqrt (2*log (2 ) ) ) ; 

energy_window  =  94:134; 
peak_height  =  500; 

time  =  20:. 05:120; 

for  k  =  1:1000; 

src  =  src_mean  +  src_std*randn ( 1 , 1 ) ; 
bk  =  bk  mean  +  bk_std*randn ( 1 , 1 ) ; 

for  i  =  1 : length (time) ; 

y  = 

(bk*time (i) )  +  (src*time (i) ) *normpdf (energy  window, energy, peak  std)  ; 

%This  it  the  histogram  of  the  photopeak  growth  through  time 

bar (energy  window, y) ; axis ( [min (energy  window)  max (energy  window)  0 

1600]); 

pause  ( . 05 ) ; 

if  max (y) -min (y)  >  peak  height; 

t ( k)  =  time ( i ) ; 
break 

end; 

end; 

end; 

t^mean  =  mean(t); 
t_std  =  sqrt (var (t) )  ; 

%  time__title  =  [ 'Mean= ',  num2str  (t_mean)  num2str  (t_std)  ]  ; 

%  hist  (t, 100 ); title (time  title); 
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Appendix  B.  Experimental  Setup  and  System  Settings 


B.l  List  of  Equipment  and  Commonly  used  settings 

Equipment  List: 

NIM  bin:  Ortec  4001 A 

Detector/PMT:  Bicron  3x3”  Nal 

High  voltage  power  supply:  Tennelec  952 A 

Pre- Amplifier:  Ortec  113 

Shaping  Amplifier:  Ortec  572 

Single  Channel  Analyzer:  Ortec  550 

Mulit-Channel  Analyzer:  Ortec  926  ADCAM 

Data  Acq.  Board:  National  Instruments  6111  E-PCI  using  a  SCB-68  interface 

Stepper  Motor:  Applied  Motion  HT23-297 

Motor  Driver:  Applied  Motion  Si  3540 

Rotatary  Table:  Velinex  B4836TS 

Motor  Controller: 

Optical  Encoder  Ring:  Renishaw  RESR  R20- 115 
Encoder  Read  Head:  RGH  20 

Decoder:  Avago  HCTL-2016  on  National  Instruments  USB  interface  6008 
Source  Positioning  System:  Velinex  Bi-slide 

Commonly  Used  NIM  settings: 

HVPS:  +800  volts 
Pre-Amp:  200  pF 
Amplifier: 

Course  Gain:  100 
Fine  Gain:  12.76 
Shaping  Time:  2  ps 

SCA: 

Mode:  Nonnal 

Lower  Window:  0.92  volts 

Upper  Window:  1.31  volts 

MCA: 

Resolution:  1024 

Source  Information 

Co-57  disk  source  with  1  cm  diameter 
Activity  Calibrated  on  20  Oct  2009  at  0.84  mCi 
AFIT  Source  ID:  T-121 
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Commonly  Used  Driver  Motor  settings: 
Motor  Resolution:  MR13 
Velocity:  VE0.1365 
Acceleration:  AC360 
Deceleration:  DE360 
Motor  Current:  CC1.5 
Step  size:  DI327 


B.2  Decoding  circuit  with  32-bit  chip  installed 


These  photos  show  the  installation  of  the  Avago  HCTC-2022  32-bit  decoder  chip  installed  to 
enable  the  processing  of  all  of  the  positional  data  without  having  to  compensate  for  the  16-bit 
chip. 
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